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The  major  objective  of  this  study  is  to  develop  a 
technique  to  analyze  delaminations  in  composite  beams  by  using 
the  finite  element  method. 

A  method  is  proposed  to  obtain  the  strain  energy  release 
rate  and  the  stresses  ahead  of  the  crack  for  the  specific  case 
of  free  edge  delaminations.  A  shear  deformable  beam  finite 
element  with  nodes  offset  to  either  top  or  bottom  side  has 
been  developed.  The  thermal  effects  originating  from  the 
curing  process  are  included  in  the  formulation  of  the  problem. 

Three  different  expressions  to  calculate  the  strain 
energy  release  rate  are  derived,  one  in  terms  of  the  force  and 
moment  resultants  in  the  elements  surrounding  the  crack-tip, 
and  the  two  others  in  terms  of  the    forces  transmitted  by  a 
rigid    element  placed  at  the  tip  of  the  crack.  The  stresses 
ahead  of  the  crack  are  obtained  by  using  discrete  spring 

ix 


elements  in  the  uncracked  portion  of  the  beam.  The  results 
calculated  are  in  good  agreement  with  results  obtained  from 
other  studies  of  the  same  problem.  The  advantage  of  using  the 
present  method  is  that  it  requires  less  computational  effort, 
and  is  capable  of  modeling  any  type  of  geometry,  loading,  and 
boundary  conditions. 

An  efficient  beam  finite  element  has  been  developed 
based  on  a  generalized  beam  theory  for  anisotropic  beams,  in 
which  the  equilibrium  equations  are  assumed  to  be  satisfied 
in  an  average  sense  over  the  width  of  the  beam.  This 
introduces  a  new  set  of  force  and  moment  resultants  for  the 
beam.  The  displacements  in  the  beam  are  expressed  in  terms  of 
three   deflections,  three  rotations  and  one  warping  term.  The 
effectiveness  of  the  method  is  illustrated  by  solving  problems 
where  the  beam  theory  solutions  are  known.  The  strain  energy 
release  rate  is  determined  for  the  cases  of  delaminated 
cantilever  beam  under  torsion  and    simply  supported 
delaminated  beam  under  transverse  loading.  Studies  are 
performed  to  understand  the  influence  of  the  position  of  the 
delamination  on  the  strain  energy  release  rate. 
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CHAPTER  1 
INTRODUCTION 

1.1  Motivation  and  Objectives 

Structural  materials  can  bear  defects  proceeding  from  the 
fabrication  processes  or  accumulate  damage  during  their 
service  life.  Frequently,  such  defects  and  damage  are  not 
visually  observable  and  have  to  be  detected  by  using  non- 
destructive evaluation  techniques  (NDE) . 

In  a  variety  of  applications,  such  as  aerospace  and 
aircraft  structures,  offshore  platforms,  bridges  and  tall 
buildings,  the  replacement  of  the  imperfect  parts  is  not  a 
cost-effective  solution  for  the  problem.  The  costs  involved  in 
keeping  the  equipment  out  of  order,  the  difficulty  of  access, 
and  other  reasons  suggest  that  there  is  a  need  for  research 
and  development  of  methods  to  predict  the  occurrence  and 
extent  of  the  damages. 

Composite  materials  and  structures,  for  example,  are 
highly  susceptible  to  damage  from  low-energy  impacts. 
Incidents  such  as  tool-drops,  improper  handling  or  hailstorms 
(aircraft  structures)  can  cause  internal  damage  and 
significant  reduction  in  compressive  strength  of  the  material. 
Under  subsequent  transverse  loadings,    the  damage  will  grow 
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causing  a  redistribution  of  the  stresses  in  the  plies  of  a 
laminate  and  consequently  reduction  in  the  material  stiffness 
and  fatigue  life.  However,  depending  on  the  subsequent  loading 
conditions  and  the  damage  characteristics,  the  structure 
often  can  still  be  used  despite  the  existence  of  damage. 

Failure  in  a  composite  material  often  initiates  in  the 
form  of  a  matrix  crack  or  a  delamination.  The  former  refers  to 
intralaminar  failure  or  transverse  crack,  while  the  latter 
refers  to  interlaminar  failure.  Both  types  of  failures  are 
illustrated  in  Fig.  l.l. 

Delamination,  or  separation  of  the  layers,  is  a  mode  of 
failure  unique  to  laminated  materials,  and  it  is  the 
predominant  failure  mode  in  continuous  fiber  composites.  Based 
on  the  location  and  direction  of  the  growth,  we  can  discern 
two  types  of  delamination:  the  free  edge  delamination  and  the 
transverse  delamination.  In  the  edge  delamination  the  growth 
initiates  at  the  load-free  edges  of  the  laminate  and 
propagates  normal  to  the  load  direction.  Under  transverse 
delamination,  also  known  as  crack  tip  or  local  delamination, 
the  growth  starts  from  a  matrix  crack  and  propagates  normal  to 
the  transverse  crack  from  which  it  originates.  The  Fig.  1.2 
shows  both  types  of  delamination.  Since  delamination  is  a  very 
commonly  observed  failure  mode  in  composite  materials,  its 
study  is  of  great  importance  in  the  field  of  failure 
prediction  for  composite  structures. 


Figure  1.1  Failures  in  composite  materials 

(a)  Matrix  crack  (b)  Delamination 


delamination 


matrix  crack 


(a) 


delamination 

(b) 


Figure  1.2  Types  of  delamination 

(a)  Free  edge  (b)  Transverse  delamination 
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1.2  Literature  Survey 

In  the  past  few  years,  several  researchers  have  studied 
delamination  occurrence  and  growth.  Usually  they  predict 
delamination  growth  from  the  strain  energy  release  rate  (G) , 
and  try  to  develop  analytical  and  numerical  methods  capable 
of  modeling  this  phenomenon.  A  brief  description  of  some  of 
these  studies  is  important  for  accurate  comprehension  of  the 
present  work. 

Chan  and  Ochoa  (1987)  developed  a  two-dimensional  finite 
element  model  to  calculate  interlaminar  stresses  and  strain 
energy  release  rate  in  composites  subjected  to  uniaxial 
tension,  bending  and  torsion.  The  element  is  an  isoparametric 
element  with  eight  nodes  and  three  degrees  of  freedom  per 
node.  The  basic  idea  in  their  work  is  to  represent  the  effect 
of  longitudinal  or  twisting  curvature  as  a  "load  effect."  The 
study  concentrates  on  establishing  relationships  of  Gj ,  G,,  and 
Gjjj  as  functions  of  crack  length  and  stacking  sequence. 

Pagano  (1974)  studied  the  case  of  fibrous  composite 
laminates  under  a  uniform  axial  extension  and  developed  an 
approximate  method  to  define  the  distribution  of  the 
interlaminar  normal  stress  along  the  central  plane  of  a 
symmetric  laminate.  A  simple  expression  for  the  interlaminar  ■ 
normal  stress  was  derived  and  the  results  obtained  by  using 
this  expression  were  compared  with  those  obtained  by  using  a 
finite  difference  approximation  of  the  elasticity  equations. 
The    two   methods    agreed   well,    and    both   results    show  that 
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significant  stresses  exist  only  in  narrow  regions  next  to  the 
free  edges.  The  length  of  this  region  is  approximately  equal 
to  the  laminate  thickness. 

Whitney  (1986)  proposed  an  approximation  to  the 
interlaminar  normal  stress  distribution  ahead  of  the  crack  and 
the  influence  of  specimen  geometry  on  Mode  I  strain  energy 
release  rate  for  edge  delamination  specimens.  The  specimens 
were  subjected  to  tension  and  compression  loadings  and  were  of 
the  classes  [0^/90^]^  or  [B/S^/e/so^]^.  Whitney  used  a  higher 
order  plate  theory,  which  considers  transverse  shear 
deformation  and  thickness-stretching.  The  effect  of  thermal 
residual  stresses  due  to  lamination  fabrication  was  considered 
in  the  formulation  of  the  problem.  The  numerical  results 
suggested  significant  influence  of  the  residual  thermal 
stresses  on  both  the  stress  distribution  ahead  of  the  crack 
and  the  Mode  I  strain  energy  release  rate.  In  general,  the 
thermal  residual  stresses  increase  the  stresses  at  crack  tip 
and  the  Mode  I  strain  energy  release  rate  for  specimens  under 
tension,  and  has  the  inverse  effect  for  specimens  under 
compression.  Again  the  results  confirmed  that  the  effect  of 
normal  stress  ahead  of  the  crack  dissipates  within  a  distance 
of  approximately  twice  the  thickness  from  the  crack  front. 

Armanios  and  Badir  (1989)  used  a  shear  deformation  theory 
to  calculate  the  interlaminar  stresses  and  energy  release  rate 
associated  with  both  Mode  I  and  mixed-mode  edge  delaminations. 
The  formulation  used  by  Armanio  and  Badir  does  not  consider 


the  thickness  strain  effect.  The  analysis  includes  the  effects 
caused  by  the  residual  thermal  stresses  and  moisture  stresses. 
The  classes  of  laminates  considered  are  the  same  ones  used  by 
Whitney  (1986) .  Results  suggest  that  the  residual  thermal 
stresses  tend  to  increase  the  interlaminar  stresses  and  total 
energy  release  rate  for  both  modes,  while  the  moisture  content 
tends  to  alleviate  the  thermal  effects.  The  value  of  moisture 
content  necessary  to  alleviate  totally  the  interlaminar 
stresses  and  the  strain  energy  release  rate  ranged  between 
0.70%  and  0.80%  in  the  examples  considered,  and  it  was  found 
that  this  value  has  a  weak  dependency  on  both  the  stacking 
sequence  and  the  mode  considered. 

Usman  and  Gandhi  (1989)  studied  the  effect  of  moisture 
and  temperature  on  the  vibrational  characteristics  of 
composite  materials.  The  nonlinear  vibration  of  moderately- 
thick  laminated  composite  plates  under  hygro-thermal 
conditions  are  considered  by  incorporating  the  effects  of 
transverse  shear  and  rotary  inertia.  The  hygro-thermal  changes 
along  the  thickness  of  the  composite  were  explicitly 
incorporated  in  the  form  of  the  deformation  field.  The 
response  of  the  composite  under  hygro-thermal  environments  was 
modeled  by  creating  one  model  formed  by  two  parts,  one 
anisotropic  solid  (composite)  and  one  ideal  fluid  (moisture) . 
The  equations  presented  in  this  work  are  significantly 
different  from  the  corresponding  equations  in  classical 
laminated  composite  plate  theory. 


Hu  et  al. (1989)  developed  a  novel  beam  element  to  study 
the  crack  propagation  in  beam  types  of  specimens.  Finite 
elements  are  used  to  model  the  specimens,  and  the  beam  element 
has  its  nodes  offset  to  one  side,  either  bottom  or  top  side. 
This  approach  allows  modeling  crack  propagation  without 
renumbering  the  nodes  as  the  crack  propagates,  and  is 
convenient  in  modeling  the  partial  contact  and  slip  at  a 
delamination  interface.  The  bottom  and  top  un-delaminated 
portions  of  the  beam  were  connected  either  by  foundation 
springs  or  rigid  elements.  For  some  specimens  the  strain 
energy  release  rate  was  computed  by  using  a  zero-volume  path 
integral  and  a  crack  closure  method,  and  the  results  were 
found  to  be  in  good  agreement.  In  the  present  study  the  beam 
theory  is  extended  to  analyze  delaminated  composite  beams. 

Sankar  (1991  b)  developed  a  beam  theory  for  laminated 
composite  beams  derived  from  the  shear  deformable  plate 
theory.  The  equilibrium  equations  are  assumed  to  be  satisfied 
in  an  average  sense  over  the  width  of  the  beam,  which 
introduces  a  new  set  of  force  and  moment  resultants  for  the 
beam.  The  displacements  in  the  beam  are  expressed  in  terms  of 
three  deflections,  three  rotations,  and  one  warping  term.  The 
closed  solutions  for  the  cases  of  a  cantilever  beam  subjected 
to  end  loads  and  a  specially  orthotropic  laminate  beam  under 
torsional  loading  are  presented.  The  solution  considers  the 
shear  deformation  and  the  edge  effect.  The  results  in  terms  of 
angle  of  twist  agree  well  with  available  solutions. 
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1.3  Scope  and  Outline  of  this  Work 

Numerical  and  analytical  methods  have  been  used  to 
express  the  complex  stress  state  near  free  edges  and  to 
calculate  the  energy  release  rate  during  crack  propagation. 
Among  several  numerical  methods,  the  finite  element  method  was 
selected  to  be  used  in  the  present  work.  Although  requiring 
intensive  computational  effort,  this  method  is  very  popular 
and  provides  accurate  solutions.  Also,  this  method  is  very 
convenient  in  modeling  structural  elements  of  variable  cross 
section. 

In  the  present  work,  the  typical  element  used  to  model 
the  structure  is  the  beam  element.  This  element  not  only  is 
very  common  in  mechanical  structures  but  also  is  the  most 
appropriate  to  model  some  standard  test  specimens. 

The  scope  of  this  work  can  be  described  as  to  develop  a 
method  to  analyze  delaminations  in  anisotropic  composite 
beams,  using  the  finite  element  method  to  model  the  problem 
and  selecting  the  strain  energy  release  rate  as  a 
delamination  parameter.  The  loading  is  considered  independent 
of  time. 

This  dissertation  consists  of  five  chapters.  Chapter  2  is 
focused  on  free  edge  delamination.  A  method  is  proposed  for 
modeling  this  type  of  delamination  using  offset  beam  finite 
elements.  The  effect  of  temperature  variation  during  the 
curing  process  (residual  stresses)  is  considered  in  the 
problem.  Chapters  3  and  4  are  devoted  to  presenting  a  new  and 


9 

efficient  method  to  analyze  composite  laminated  beams,  using 
a  more  complete  and  accurate  beam  finite  element.  The 
capability  of  this  new  finite  element  in  modeling  the  problem 
was  verified  by  solving  the  problem  of  a  cantilever  beam 
subjected  to  different  end  loading  conditions,  including 
torsion.  Applications  to  delaminated  beams  under  torsion 
loading  and  the  simulation  of  a  quasi-static  impact  test  on  a 
delaminated  beam  are  also  presented  in  these  chapters. 
Finally,  chapter  5  summarizes  the  principal  conclusions  of  the 
present  study,  and  presents  suggestions  for  future  studies  in 
this  field. 


CHAPTER  2 

AN  OFFSET  BEAM  FINITE  ELEMENT  FOR  MODELING 
FREE  EDGE  DELAMINATIONS  IN  COMPOSITES 

2 . 1  Introduction 

The  separation  of  plies  or  delamination  is  a  predominant 
failure  mode  in  laminated  composite  structures.  It  occurs  when 
interlaminar  stresses  exceed  the  interlaminar  strength  of  the 
material,  or  as  a  result  of  the  fabrication  process, 
environmental  degradation  and  improper  handling  of  the 
material.  Common  examples  of  such  failure  regions  during  the 
normal  use  of  composite  structures  are  those  near  cut-outs, 
holes  and  free  edges.  In  these  regions,  the  abrupt  changes  in 
the  geometry  or  material  properties  tend  to  intensify  the 
stresses  when  the  laminate  is  subjected  to  some  external 
loads. 

Specifically,  the  occurrence  of  delamination  in  free 
edges  has  been  receiving  increasing  attention  of  many 
investigators,  in  their  efforts  to  understand  and  to  prevent 
delaminations  in  composite  structural  designs.  The  reason  is 
the  presence  of  high  interlaminar  stresses,  the  "peeling" 
stress,  in  the  neighborhood  of  a  free  boundary.  The  "free  edge 
problem"  usually  is  modeled  by  a  plate  subjected  to  a  uniform 
axial  strain  (  e^  )  ,  as  shown  in  Fig.  2.1. 
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(applied  strain) 


Figure  2.1  Free  edge  delamination 
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During  the  usual  fabrication  process,  the  composite 
material  is  subjected  to  temperatures  of  several  hundred 
degrees  Fahrenheit,  and  then  it  is  cooled  down  to  the  room 
temperature.  The  composite  is  assumed  to  be  stress-free  at  the 
curing  temperature.  The  change  in  temperature  causes 
unavoidable  stresses  to  the  composite  material  layers  when 
cooled  to  room  temperature.  Usually  these  stresses  are  called 
residual  stresses  or  thermal  stresses,  and  they  are  considered 
in  the  formulation  of  the  present  problem. 

A  simple  procedure  is  developed  in  this  study  to  obtain 
an  approximate  solution  for  the  peel  stresses  ahead  of  the 
crack  tip  and  to  calculate  the  strain  energy  release  rate,  for 
the  case  of  free  edge  delaminations  including  thermal  effects. 
The  problem  depicted  in  Figure  2.1  is  modeled  using  a  shear 
deformable  beam  finite  element  with  offset  nodes  (Sankar, 
1989) .  The  three-dimensional  problem  is  modeled  as  a  plane 
problem  in  the  x-z  plane  assuming  that  the  stresses  and 
strains  do  not  vary  along  the  y-direction.  This  will  be  true 
for  most  of  the  length  of  the  plate  in  the  y-direction.  The 
delamination  plane  divides  the  beam  into  two  parts,  the  top 
sublaminate  and  the  bottom  sublaminate  (Figure  2.2),  with 
nodes  offset  to  the  bottom  and  top,  respectively.  The 
delamination  plane  is  assumed  to  be  the  weakest  one,  therefore 
the  crack  propagates  along  this  plane,  without  crack 
branchings . 


bottom  sublaminate 


Figure  2.2  Sublaminates 
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When  calculating  the  energy  release  rate,  the  uncracked 
region  of  the  beam  is  connected  by  common  nodes,  except  the 
tip  of  the  crack,  where  a  rigid  element  is  used  to  connect  the 
top  and  bottom  nodes  (Fig.  2.3).  This  rigid  element  not  only 
ensures  the  continuity  of  displacements  and  rotations  at  that 
section,  but  also  brings  in  its  formulation  the  generalized 
forces  that  are  transmitted  between  the  sublaminates .  These 
forces  are  used  in  the  calculation  of  the  energy  release  rate. 
In  evaluating  the  peel  stresses  ahead  of  the  crack,  the 
uncracked  region  of  the  beam  is  connected  by  spring  elements 
(Fig.  2.4) ,  stiffness  of  which  has  to  be  judiciously  selected. 

The   strain   energy  release   rate   has   been  successfully 
adopted  by  researchers  in  characterizing  delamination.  In  this 
study   its   value    is   obtained   from   three   different  methods 
(Sankar    and    Pinheiro,    1990)  .    The    first    one    uses    in  its 
formulation   the   forces   transmitted   by  the   crack-tip  rigid 
element.  In  the  second  method,  the  result  is  obtained  by  using 
a  method   similar  to  the  crack-closure  method  used   in  two- 
dimensional  fracture  problems.   Finally,   composites  are  very 
brittle  materials, so  that,  the  plastic  zone  in  the  vicinity  of 
the  crack  tip   is  negligible,    and  hence  the   linear  elastic 
fracture  mechanics  can  be  applied  in  which  case  the  J-integral 
is  equivalent  to  the  energy  release  rate.  Therefore,  the  third 
method  consists   in  obtaining  the  J-integral   from  the  beam 
model.     The    force    and    moment    resultants    in    the  region 
surrounding  the  crack-tip  are  used  in  this  calculation. 
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cracked  region     uncracked  region 


Figure  2.3  Finite  element  model  for  strain  energy 
release  rate 


Figure  2.4  Finite  element  model  for  peel  stresses 
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The  results  from  the  three  different  methods  are  compared  and 
found  to  be  in  good  agreement. 

2 . 2  Laminated  Beam 

In  this  section  a  homogeneous  orthotropic  laminate 
subjected  to  a  uniform  initial  strain  =  e^  applied  over  a 
length  2L  is  analyzed.  The  thickness  of  the  laminate  is  2h  and 
the  width  is  2b.  The  thickness  2h  is  much  smaller  than  2L  and 
2L  much  smaller  than  2b  (Fig.  2.1).  Due  to  symmetry 
conditions,  just  one  half  of  the  laminate  is  considered  here. 

The  problem  is  divided  into  two  different  phases.  In  the 
first  phase  no  initial  strain  is  applied  to  the  specimen, 
which  is  subjected  only  to  thermal  stresses  due  to  the  curing 
process.  In  the  second  phase  the  cured  specimen  containing  the 
edge  delamination  is  subjected  to  the  axial  strain  =  in 
addition  to  the  residual  stresses. 

The  equations  derived  in  the  next  sections  refer  to  a 
laminated  beam  situated   above  the  delamination  plane.  The 
expressions    for    the     laminated    beam    situated    below  the 
delamination  plane  differ  only  in  the  integration  limits 
(  -h  to  0 ,  instead  of  0  to  h  ) . 

2.3  Fabrication  Process 

During  the  composite  fabrication  process  there  are 
heating  and  cooling  phases.  In  the  cooling  phase,  the 
temperature  drops  from  the  stress-free  temperature  to  the  room 


17 

temperature.  In  this  case,  each  layer  influences  the  expansion 
or  contraction  of  the  other  because  of  the  difference  in  their 
coefficient  of  thermal  expansion.  A  similar  situation  exists 
when  the  composite  absorbs  moisture  during  service.  As  a 
result  of  the  temperature  variation,  each  layer  has  residual 
stresses  after  the  laminate  fabrication. 

The  effects  of  the  residual  stresses  on  the  laminate  as 
a  whole  are  calculated  and  expressed  in  terms  of  thermal 
forces  and  moments.  In  the  next  phase,  these  thermal  forces 
and  moments  are  added  to  the  external  forces  and  moments. 

During  the  cooling  process  the  in-plane  and  transverse 
displacements  are  assumed  to  be  of  the  form: 

u{x,z)-  Uq{x)  +       j^ix)  (2.1) 
v{y,z)-  v^iy)  +  z^  yiy) 
w(x,  z)  -  Wq  ix) 

where  Ug  and  Vp  are  the  axial  displacements  of  points  on  the 
x-y  plane  in  the  x  and  y  directions  respectively,  and  and 
are  the  rotations  about  the  y  and  x  axes  respectively.  It 
is  assumed  that  the  modification  in  the  thickness  of  the 
laminate  is  negligible  when  compared  with  the  modifications  in 
the  other  two  dimensions.  Therefore,  w  is  assumed  to  be 
constant  through  the  thickness  of  the  laminate  (no  thickness 
strain) . 
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From  Equations  2.1  the  strains  for  the  cooling  process 
are  derived  as 


(2.2) 


A  state  of  plane  stress  normal  to  the  z  direction  is 
assumed  during  the  curing  process.  Therefore,  the  thermal 
stresses  developed  within  the  k*^  layer  (Agarwal  and  Broutman, 
1980)  are 


'^11^12^16 
^12^22^26 
A6^26^66 


XX      »  X 

e  -C 
yy   '  y 


-c 


xyj 


(2.3) 


Where  Q..  are  the  stiffness  coefficients  and  can  be  found  in 


Whitney  (1987),  and 


Cx-  <^x^T 


(2.4) 


In  the  above  equations  a  is  the  thermal  expansion 
coefficient  and  AT  is  the  change  in  temperature. 

By  definition,  the  force  (N)  and  moment  (M)  acting  on  the 
top  laminate  are 


.,M^)  -Jo ^[1,  z]  dz 


0 

h 


{Ny,My)       -fay[l,Z]  dZ 


(2.5) 
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In  the  present  case,  there  are  no  net  resultants  acting 
on  the  laminate.  Therefore,  considering  the    definition  of 
force  and  moment  resultants,  we    have  for  the  whole  laminate 


dz 


(2.6) 


and 


MA 


i-s/C;), 


zdz 


(2.7) 


where  n  is  the  total  number  of  layers  on  the  top  sublaminate 
and  (h,^  -  h,^.^)   is  equal  to  the  thickness  of  the  k**'  layer 
(Fig.  2.5). 

layer  number 


T 


n 


n-1 


plane  of  delamination 


Figure  2.5  Values  of  h  for  top  sublaminate 
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The   stiffness  coefficients  of  the  top  sublaminate  are 
defined  as 


(2.8) 


i/D  ~  1/2,. ..,6 

In  the  above  equation,  a  shear  correction  factor  equal  to 
5/6  (Whitney,  1987)  is  introduced  in  the  calculation  of  the 
terms  ajj,b..  and  d,.,.,  wherever  the  index  i  is  greater  then  3. 

Substituting  Eq.  2.3  in  Eq.  2.6  and  2.7,  and  using  the 
definition  of  the  stiffness  coefficients  in  Eq.  2.8,  we  obtain 
the  following  expressions  for  the  force  and  moment  resultants 


xf 


-  (o°)  - 


•12 


^la 


V^12  C?i2j 


/_0 


(2.9) 


^12  ^12 


^22  ^22 
^22/ 


■y  - 


(2.10) 


where  N^^  and  M^^  are  the  thermal  force  and  moment  in  the  x 
direction  and  so  are  N^^  and  M^^    in  the  y  direction. 
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The  thermal  forces  and  moments  are  defined  as 


[nJ,mJ]-  fiO^^C^^  Q^^Cy^QrsC^cy)  n,z]dz  (2.11) 


[<,m/]-  /(Cx2C,  +  022Cy  +  02sC^)  dz  (2.12) 

0 

Defining  the  matrices 


and 


[L]- 


[M]' 


^11  ^11^ 


\,'^22  "^22^ 


Equations  2.9  and  2.10  take  the  form  of 


and 


/0\ 


)  -  [L] 

+  [M] 

)  -  m 

From  Eq.  2.17  we  obtain 


-  [N] 


(2.13) 


(2.14) 


(2.15) 


(2.16) 


(2.17) 


(2. 18) 


and,  after  substituting  from  Eq.  2.18  in  Eq.  2.16,  we  obtain 

(2.19) 


■  T 


[p]  r +  [M]  m  -1 


/  7^ 
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where 

[P]-  [L]  -  m  m-"-  IM] 

Eq.  2.19    defines  the  thermal  force  and  moment,  and 
respectively,  as 

-  [M]  m-^l^'[P]{\A  (2.20) 

From  now  on,  the  thermal  force  and  moment  resultants 
defined  in  Eq.  2.20  are  considered  as  external  forces  and 
moments  acting  on  the  composite  specimen  (Agarwal  and 
Broutman,  1980) .  The  finite  element  method  is  used  to  solve 
Eq.  2.20  and  to  obtain  the  solution  for  e°^(x)  and  K^ix)  .  For 
bottom  element  the  derivation  is  the  same,  except  for  the 
limits  of  integration  in  Eq.  2.5,  2.8,  2.11  and  2.12. 

2.4  Free  Edge  Problem 

In  defining  the  in-plane  and  transverse  displacements  for 
the  free  edge  problem  (second  phase),  Eq.  2.1  remains  the 
same,  except  for  displacement  in  the  y  direction,  which  takes 
the  form  of 

v(y,z)  -  e^y  (2.21) 
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Where  e^,  represents  the  applied  strain  in  the  y  direction. 
Consequently,  the  strain  e^^  in  Eq.  2.2  is  replaced  by 

e^y  -  (2.22) 

Because  the  dimension  2b  in  Figure  2 . 1  is  considered  much 
greater  than  the  dimension  2L,  it  is  assumed  that  a  state  of 
plane  strain  parallel  to  the  x-z  plane  exists  in  the  second 
phase  of  the  free  edge  problem.  The  stress-strain  relations 
for  the  k'^*'  layer  in  this  case  are 


0  0 


0 

(e,  y 

XX 

0 

y  xz, 

(2.23) 


The   in-plane  force  and  moment  acting  on  the  beam  are 
shown  in  Fig.  2.6 


Figure  2.6  Force  and  moment  resultants  (top  sublaminate) 
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P,  M  and  V  are  defined  as 

h 

[P,M,V]-  [  [a^,ZO^,X^^]   dz  (2.24) 

0 

or,  using  the  definitions  given  in  Eq.  2.8  and  2.23,  as 

(2.25) 


^11 

0  ^ 

XX 

\m 
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1° 

0 

y  xz, 

[o) 

Where  and  are  the  force  and  moment  induced  by  the 
applied  strain  e^.  They  are  defined  as 

h 

[P^^^^]- 1^(1,2)  dz  -  [ai2,23i2leo  (2.26) 

0 

The  thermal  force  and  moment  resultants  defined  in  Eq. 
2.20  are  added  to  the  problem  in  the  form  of  external  forces, 
and  hence  the  problem  takes  the  final  form  of 


(2.27) 
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(2.28) 


The  last  equation  can  be  expressed  in  a  more  concise  form 


as 


F  -  S  e  -  R  (2.29) 
where  S  is  the  stiffness  matrix,  components  of  which  are 


[5]- 


b 

b. 
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11  c?ii  0 
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(2.30) 


and  R  is  the  vector  that  represents  the  forces  and  moments 
induced  by  the  applied  strain  and  by  the  thermal  effects 


{Ry 


pH 

[o  ) 

lo  j 

(2.31) 


The  vector  F  refers  to  the  external  forces  acting  on  the 


beam  and  its  components  are 


(F)  - 


M 


(2.32) 
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while  e  is  the  vector  of  deformations  in  the  form  of 


(2.33) 


Finally,  the  Equation  2.27  can  also  be  expressed  as 


e  -  c (F+R) 


(2.34) 


where  c  =  s'^    and  is  the  compliance  matrix. 

2.5  Finite  Element  Solution 

The  Finite  Element  Method  is  used  to  solve  Eq.  2.29,  and 
a  beam  finite  element  with  offset  nodes  is  selected  to  model 
the  problem.  The  derivation  of  this  specific  beam  finite 
element  was  done  by  Sankar  (1989).  The  finite  element  has  two 
nodes  and  three  degrees  of  freedom  (\x,(f/  and  w)  at  each  node. 
The  nodal  forces  and  moment  for  the  i***  node  are  denoted  by  p,., 
m,.  and  v,. .  In  the  present  study  a  slightly  different  method  is 
used  to  derive  the  stiffness  matrix. 

The  average  force  and  moment  resultants  (P,  M  and  V)  at 
the  center  of  the  element  are  given  by 


(P\  ((P2  -  Pi)  /  2\ 
M  -  m^)  /  2 

W)    [iv^  -  V,)  /  2^ 


(2.35) 


The    deformations    at    the    center    of    the    element  are 
expressed  in  terras  of  the  nodal  displacements  as 


(e)- 


zx) 


(^(T^  +  T,)  /2+  IL] 


(2.36) 


where  L  is  the  element  length. 

The  beam  finite  eleraent  forces  and  degrees  of  freedom  are 
shown  in  Fig. 2. 7.  By  substituting  from  Eq.  2.33-2.34  into  Eq. 
2.29,  the  relations  between  nodal  forces  and  displacements  are 
obtained  in  the  form  of 

f  -  k  g  -  R  (2.37) 
where  the  vectors  of  nodal  forces  (f)  and  displacements  (q) 

have  the  following  components 


(g)  ^-  [Ui,l|fi,  W^,  U2,^2'^2^ 


(2.38) 


The  vector  R  and  its  components  are  defined  by  Eq.  2.20, 
2.26  and  2.31.  The  elements  of  the  stiffness  matrix  k  are 
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2 . 6  The  Finite  Element  Models 

Fig.  2.3  and  2.4  show  the  models  used  in  modeling  the 
free  edge  delaminations .  They  consist  of  two  sublaminates 
above  and  below  the  delamination  plane.  The  assumptions  are 
that  the  delamination  plane  is  the  weakest  plane  in  the 
material  and  that  there  is  no  frictional  resistance  against 
sliding  of  the  sublaminates.  The  first  assumption  means  that 
the  crack  will  propagate  along  the  delamination  plane,  without 
crack  branching. 

In  the  uncracked  region,  the  top  and  bottom  portions  of 
the  beam  are  connected  by  common  nodes,  when  calculating  the 
strain  energy  release  rate  (Fig.  2.3)  ,  and  by  spring  elements, 
when  predicting  the  stress  ahead  of  the  crack  (Fig.  2.4). 
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(a) 


Figure  2.7  Force  and  degrees  of  freedom 

(a)  Nodal  and  average  forces  (b)  Degrees  of  freedom 
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2.7  Strain  Energy  Release  Rate 

In  the  calculation  of  the  strain  energy  release  rate  (G)  , 
one  rigid  element  is  placed  in  the  tip  of  the  crack  connecting 
both  sublaminates  (Fig.  2.3).  The  rigid  element  has  nine 
degrees  of  freedom.  Apart  from  the  three  displacements  at  each 
node,  the  three  generalized  forces,  F^,  transmitted  by  the 
rigid  element  are  also  considered  as  unknown  degrees  of 
freedom.  The  constraint  equations  corresponding  to  the  rigid 
element  are 


'O3     O3  -J3^ 

03  03  J3 


(2.40) 


where  F,  and       are  the  nodal  forces  at  nodes  1  and  2  and 
and        are  the  nodal  displacements.   The  symbol  O3  represents 
the  3x3  null  matrix  and  I3  the  3x3  identity  matrix. 

Considering  a  zero  volume  path  surrounding  the  crack  tip 
as  shown  in  Fig.  2.8,  the  expression  for  the  J-integral 
relative  to  two-dimensional  states  can  be  found  in  Hellan 
(1984)  as 


/du  ■ 
iwdx^-Pi^ds) 


(2.41) 
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where 

Pi  -  rijOji  i-1,2 

For  the  present  case,  and  considering  the  path  r  very 
thin  (zero-volume),  the  Eq.2.41  takes  the  form  of 

J  -  fwdz  -  fn^o^u^ds-fn^x^^w,,ds  (2.42) 
r  r  r 

where  W  is  the  strain  energy  density  and  is  defined  as 

-|  (Oxx^xx  +  Oyyfiyy       ^z;cYzx)  (2.43) 

The  path  r  has  four  segments,   one  in  each  sublaminate 
behind  and  ahead  of  the  crack  tip  (Fig.  2.8).  The  values  of  n^^ 
(unit  vector  normal  to  the  path)  and  ds  for  each  element  are: 
Elements  1  and  4:  n^^  =-1;  ds  =-dz 

Elements  2  and  3 :  n^^  =  1 ;  ds  =  dz 

Considering  the  Element  1  in  Fig. 2. 8  we  have: 

h  h  h 

J^-  -fwdz^  ja^u,,dz^jx^^w_^dz  (2.44) 

0  o  0 

where  the  subscript  stands  for  the  element  number. 

The  first  term  of  the  right  hand  side  of  Eq.2.44  is 
defined  by  using  Eq.2.43  and  the  stress-strain  relations  given 
in  Eq.2.2  3  as 

ti  h  (2.45) 

0  0 
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Figure  2.8  Zero-volume  path  for  J-integral 
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From  the  definition  of  the  stiffness  coefficients  and  Eq. 
2.30  and  2.3  3  we  obtain  the  following  relation 

h 

e^Se-  f  {Ai^^'^Qs^yz.')  dz  (2.46) 

0 

Therefore,  the  Eq.2.45  can  be  modified  and  expressed  as 

h 

jwdz-^e^Se^  +  a^^e^G^  ^  ^a,,e,'  (2.47) 

The  second  term  of  the  right  hand  side  of  Eq.2.44  is 
modified  in  a  similar  way  to  what  has  been  done  for  the  first 
term. 

h  h 

{^xxU.x<^z-f  (^iifixx^  +  ^iaG^eJ  dz  (2.48) 

0  0 

or 

h  h 

f^^^,:c<^z-f  (^iie«'-^^i2G«.eo  +  ^55Yzx'-^55Yz/)  dz  (2.49) 

0  0 

Finally, 

h 

I'^xxU.xdz-e^^Se^  +  a^^e.e^-a^^yJ  (2.50) 
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The  third  term  of  Eq.  2.42  is 

h 

f-^zxK.^Z-^x^i  (2.51) 

0 

where       and  ip^  are  the  shear  force  resultant  (Eq.  2.24)  and 
the  rotation  at  the  crack  tip,  respectively,  for  the  element 
number  1. 

Substituting  Eq.   2.48,   2.50  and  2.51  into  Eq.    2.44  we 
obtain 

j^-le^-'Se^  +  (2.52) 
where  U^.  is  a  energy  term  defined  by 

-  -|  322^02-355  Yz;c'  (2.53) 

similar  expressions  can  be  obtained  for  the  other 
elements  in  the  form  of 

-  -\e^^Se^  -  (2.54) 

For  linear  elastic  fracture,  the  J-integral  is  equivalent 
to  the  energy  release  rate  (Hellan,  1984) .  This  is  valid  for 
fractures  in  brittle  composite  materials,  where  the  plastic 
zone  in  the  vicinity  of  the  crack  is  negligible.  Hence, 
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G 


-  J  -       +  J2  +  J^i  +  - 


(2.55) 


The  Eq.2.55  can  be  calculated  in  a  more  convenient  form, 


where  the  forces  F.  (1=1,2,3,4)  are  calculated  from  the  finite 
element  solution,  and  so  are  the  compliance  matrices  c^  and  c^^ 
from  the  input  data. 

The  elements  2  and  3  are  in  the  uncracked  region,  and  so, 
they  are  under  the  same  deformations.  We  thus  obtain  the 
relations, 

cP-cF  (2.57) 

Substituting  the  last  equation,  and  also  the  equilibrium 
relations  (F^  +  F^  =  Fj  +  F3)  ,   in  Eq.  2.56  we  obtain 


G-^  (F/c,F,  +  fJc^F,  -  F/c.F^  -  F.^C^F,) 


(2.56) 


G-1  {F,^c,  -  F/c^)  (F.-FJ 


(2.59) 


or,  using  again  the  equilibrium  relations. 


(2.60) 


36 

Actually,   (F,  -  Fj)  is  the  force  transmitted  by  the  rigid 
element,  F^,  between  bottom  and  top  crack  tip  nodes. 
Therefore,  the  strain  energy  release  rate  (G)  can  be  expressed 
in  a  more  convenient  form  as 

G  -  i.F/(c,  +  c^)  (2.61) 

where  F^  are  considered  unknown  degrees  of  freedom  and  part  of 
the  finite  element  solution  of  the  problem  (Eq.  2.40). 

One  can  also   evaluate  the  value  of  the   strain  energy 
release  rate  using  a  method  analogous  to  the  virtual  crack 
closure  method  in  two-dimensional  fracture  problems 
(   Hellan,    1984) .    In  this  method,    the  energy   necessary  to 
"close"  the  crack  is  equal  to  the  energy  release  rate. 

In  Fig.  2.9,  as  the  crack  propagates,  the  crack  tip  Nodes 
1  and  2  have  displacements  and  Ag,  respectively,  and  the 
crack  tip  moves  to  Node  5.  If  the  length  of  the  elements 
behind  the  crack  tip  is  sufficiently  small,  AL,  then  the  crack 
opening,  (Aq^-Aqj)  ,  will  be  approximately  equal  to  the 
difference  in  displacements  of  Nodes  3  and  4,(Aq3-Aq^). 
Therefore,  the  problem  is  analyzed  considering  the 
displacements  of  Nodes  1  and  2  as  the  tip  displacements  of 
cantilever  beams  of  length  AL  under  applied  forces  F^  on  each 
beam.   For  very  short  cantilever  beams   (thickness/ length  «1) 


Figure  2.9  Crack-closure  method 
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the  tip  deflections  (q^jp)  can  be  calculated  using  the 
expression 

g^ip  -  L  [c]  F  (2.62) 

where  L  and  F  are  the  beam  length  and  the  forces  applied  at 
the  tip  of  the  beam,  respectively,  and  c  is  the  compliance 
matrix.  In  the  present  case, 

Ag,--ALc,F,  <2.63) 


or 


Ag^-Ag^ 


-(c,  +  c^)F^  (2.64) 


AL 

Substituting  the  last  equation  in  Eq.  2.61  gives 

^-  lil^r^(Ag,-Ag,)  (2.65) 

In  this  equation,  the  displacements  in  the  tip  nodes  and  the 
forces  in  the  rigid  element  are  part  of  the  finite  element 
solution. 

In  the  present  study,  the  strain  energy  release  rate  was 
calculated  by  three  different  methods.  These  methods  are 
represented  by  Eq.2.56,  2.61  and  2.64,  and  the  results 
obtained  from  the  three  equations  are  in  very  good  agreement. 


39 

2.8  Stresses  Ahead  of  the  Crack 

In  evaluating  the  stresses  ahead  of  the  crack  tip,  the 
uncracked  region  was  connected  by  discrete  spring  elements 
(Fig.  2.4).  The  stress  ahead  of  the  crack  in  the  i^^  element 
was  obtained  by 

where  V  and  L  are  the  vertical  force  in  the  spring  element  and 
the  length  of  the  element,  respectively.  The  vertical  forces 
are  obtained  from  the  Finite  Element  solution  of  the  problem, 
and  the  element  lengths  are  part  of  the  input  data. 

The  determination  of  the  spring  constant  for  the  spring 
elements  is  of  significant  importance  in  appropriately 
modeling  the  problem.  The  spring  constants  are  derived  from 
the  change  in  thickness  of  a  beam  subjected  to  transverse 
forces.  The  exact  change  in  thickness  is  obtained  from  the  two 
dimensional  elasticity  theory  (Timoshenko,  1970)  .  In  the  first 
step,  just  the  top  part  of  the  beam  is  considered  and  it  is 
represented  in  Fig.  2.10.  From  the  elasticity  theory  we 
have,   for  this  case 

"T^  (y  y^-c^y-|c3)  (2.67) 


Figure  2.10  Model  for  spring  constant  calculation 
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where  E*  is  defined  as 


^"33-  (2.68) 


and 


A-     ^-^12V21-V23V32-V3lVl3-2V2iV32Vl3  ^^'^^^ 

^1  ^2  ^3 

The  vertical  displacement  (v)  is  obtained  by  integrating 
the  deformation        in  the  y  direction.  The  result  is 


where       is  a  constant. 

From  Eq.  2.70  the  displacements  at  the  positions  y=-c  and 
y=0  are  determined  as 

v(-c)-  ^  Cl  (2.71) 

16 

v(0)-  CI 

The  relative  displacement  between  these  two  particular 
positions  is 

V,-  W-O-WO).  -M^  (2.72) 
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For  a  beam  with  unit  area     we  have  that  P  =  q.   In  the 
present  case,     c  is  equal  to  half  the  thickness  of  the  beam 
(c  =  h/2).  Therefore  Eq.  2.72  is  rewritten  as 


where  k  is  the  spring  constant. 

By  comparing  between  Eqs.  2.73  and  2.74  we  have  the  value 
for  the  spring  constant  per  unit  length  considering  just  the 
top  part  of  the  element. 

k  -  IMl  (2.75) 

13h  ^  ' 

Considering  both  parts,    top   and   bottom,    we   have  two 
springs  placed  in  a  series  arrangement  (Fig.   2.10).  The 
spring  constant  equivalent  (k^^)  to  this  kind  of  arrangement 
is  given  by 


P  - 


13h 


(2.73) 


For  a  spring  model  we  have 


P-  kV^ 


(2.74) 


k. 


eg 


(2.76) 
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where  the  subscripts  t  and  b  stand  for  top  and  bottom 
respectively.  Substituting  Eq.  2.80  in  the  last  equation  we 
obtain  the  final  expression  for  the  spring  constant  used  in 
the  present  study 


2.9  Numerical  Results 

Some  results  from  the  investigation  of  strain  energy 
release  rate  and  the  interlaminar  stresses  in  a  free  edge 
specimen  can  be  found  in  Whitney  (1987) ,  (Pagano,  1989)  and 
Raju  (1986) .  These  researchers  computed  results  for  the  same 
composite  specimen.  The  present  study  uses  the  same  specimen 
to  make  the  comparison  easier. 

The  free-edge  specimen  width  to  thickness  ratio  (b/h)  is 
25  (Fig  2.11-2.16),  which  is  typical  for  conventional  edge 
delamination  test  (EDT)  specimens.  Considering  a  coordinate 
system  where  and  X2  are  parallel  and  transverse  to  the 
fibers  respectively,  and  Xj  is  through  the  thickness,  the  ply 
properties  have  the  following  relations  and  values: 


32E* 


(2.77) 


13  (ht  +  Afc) 


E1/E2     =  14 


E2/E3     =  1 


G^j/Ej  =  0.533 


G23/E2  =  0.323 


=  0.3 


=  0.55 


a 


=  -9  E-07/°C 


Otg  =  "3  =  23  E-06/°C 
yffj  =  yffj  =  5560  E-06/% 


0:  =  0 


'  s 


'  s 
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where  E.  is  the  elastic  moduli  in  the  i*^  direction,  G..  is  the 
shear  modulus  in  the  x,.-Xj.  plane  and  v..  are  the  Poisson's 
ratios.  Finally,  the  thermal  expansion  coefficients  and  the 
moisture  swelling  coefficient  in  the  x,.  direction  are  denoted 
by  a.  and  0.^  respectively.  These  properties  and  relations 
among  them  are  typical  of  current  high  performance 
graphite/ epoxy  unidirectional  composites. 

The  specimens  considered  in  the  present  study  have  the 
following  lay-ups: 

Specimen  A:   [903/03]^  Specimen  B:  [-eo/eOgZ-eo/Oj], 

Specimen  C:   [-45/452/-45/02]g      Specimen  D:  [O2/-6O/6O2/-6O], 

A  reasonable  stress-free  temperature  for  standard 
graphite/ epoxy  laminates  is  138°C  (280°F)  (Whitney,  1987) . 
Assuming  21°C  (70°F)  as  the  room  temperature,  the  difference 
in  temperature  (AT)  that  causes  residual  stresses  in  the  plies 
is  equal  to  -117°C  (210°F) . 

The  critical  axial  strains  (eg)  experimentally  determined 
in  edge  delamination  tests  are  0.3%  for  angle-ply  laminates, 
and  1%  for  the  special  case  of  a  bidirectional  class  of 
laminates  (Whitney,  1986) .  Therefore,  these  values  of  strains 
are  used  in  the  numerical  examples. 

2.9.1  Strain  energy  release  rate 

The  normalized  strain  energy  release  rate,  G^,  is  given 

by 


In  the  present  study  the  energy  release  rate  was 
calculated  using  Eqs.  2.56,  2.61  and  2.65.  The  results 
obtained  from  these  three  methods  were  ^n  good  agreement  and 
they  are  refered  to  as  "Present"  in  Table  2.1.  Comparison  of 
the  results  from  the  present  study  with  available  results  is 
also  present  in  Table  2.1.  The  relative  differences  between 
the  results  from  the  offset-node  method  and  each  of  the  other 
three  approaches  are  shown  in  Table  2.2. 

From  Tables  2.1  and  2.2  it  may  seen  that  the  present 
results  compare  very  well  with  Raju  (1986) ,  whereas  the 
comparison  with  the  other  two  methods  is  not  good  in  a  few 
cases.  It  should  be  mentioned  that  of  the  three  methods  we 
compared  with,  Raju  (1986)  used  a  detailed  finite  element 
analysis. 

2.9.2  Stress  ahead  of  the  crack 

The  interlaminar  stresses  ahead  of  the  crack  tip  are 
evaluated  using  spring  elements  in  the  uncracked  region  of  the 
specimen.  The  spring  constant  for  these  elements  is  calculated 
by  Eq.  2.77  and  found  to  be  equal  to  300x  lo'  N/m  . 
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Fig.    2.11-2.15  show  the  results  for  peel  stresses  for 
different  laminate  configurations,   with  and  without  thermal 
effects.  From  the  results  we  find  that  the  effect  of  residual 
thermal  stresses   is  to  cause  increase   in  the  maximum  peel 
stress  at  the  free  edge.  The  effect  is  significant  for 
specimen  A  (   [9O3/O3]   )  and  less  significant  for  specimen  C 
(  [-45/452/-45/0213  ).  In  Fig.  2.11  the  stress  distribution  is 
compared  with  Whitney  (1987)  for  A=  0°C  and  A=  -117 °C  for  the 
specimen  B  (  [-eo/eOjZ-eo/Ojlg  ).  The  agreement  is  excellent. 
For  the  bidirectional  laminates,  specimen  A,  the  influence  of 
temperature  can  be  better  appreciated  in  Fig. 2. 15,  where  the 
temperature  variation  goes  from  zero  to  -120 °C  in  increments 
of  -30°C. 

Besides  thermal  stresses,  some  studies  on  effect  of 
moisture  absorption  on  the  free  edge  stresses  were  also 
performed.  In  fact,  the  analysis  procedure  is  similar  except 
that  a  AT  has  to  be  replaced  by  )S  AH,  where  yff  is  the  swelling 
coefficient  and  AH  represents  the  variation  in  the  moisture 
concentration.  Fig.  2.16  shows  the  influence  of  the  moisture 
level  variation  on  specimen  A.  The  moisture  level  was  allowed 
to  vary  from  0%  to  1.2%  of  the  laminate  weight,  which  reflects 
feasible  conditions.  This  result  agrees  with  those  obtained  by 
Armanios  and  Badir  (1989).  Unlike  the  thermal  stresses,  the 
residual  moisture  tends  to  decrease  the  interlaminar  stress 
close  to  the  crack  front. 
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2.9.3  Conclusions 

The  residual  thermal  stress  tends  to  increase  the  total 
energy  release  rate  and  the  interlaminar  stress  at  the  tip  of 
the  crack,  while  the  residual  moisture  alleviates  these 
effects. 

The  results  obtained  from  the  offset  nodes  method  are 
very  close  to  those  obtained  by  Whitney  (1987),  Pagano ( 1989) , 
Raju(1986)  and  Armanios  and  Badir  (1989)  using  different 
methods.  Therefore,  the  present  offset-node  method  is  a  good 
method  to  predict  the  energy  release  rate  and  stresses  ahead 
of  the  crack  tip  in  the  free-edge  problems.  The  present  finite 
element  method  has  the  facility  to  model  any  type  of  geometry, 
loading  and  boundary  conditions,  and  avoids  renumbering  the 
nodes  and  elements  as  the  crack  propagates. 

The  offset  node  method  can  be  used  in  studying  dynamic 
delamination,  and  can  be  extended  to  problems  of  sublaminate 
buckling  and  delaminations  in  composite  plates. 
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CHAPTER  3 

MODELING  DELAMINATION  IN  COMPOSITE  BEAMS 
AND  APPLICATION  TO  TORSION  PROBLEM 

3 . 1  Introduction 

The  application  of  composite  materials  is  increasing  day 
by  day.  Not  only  traditional  areas  like  aerospace  and 
automobile  structures,  but  also  several  other  areas,  such  as 
robotics,  marine  industries  and  medical  devices  and 
prostheses,  have  shown  that  the  application  of  composite 
materials  is  major  factor  in  the  improvement  of  their 
mechanical  structures. 

The  advances  in  the  manufacturing  processes  of  composite 
materials,  the  development  of  new  matrices  and  fibers,  and  the 
progress  in  methods  for  study  of  the  mechanical  behavior  of 
composites  have  been  made  it  possible  for  the  composite 
materials  to  replace  traditional  materials  in  several 
mechanical  structures.  Usually  composite  materials  are 
considered  as  a  possible  and  adequate  solution  whenever  the 
weight  savings  and  the  capability  of  "tailoring"  the  material 
in  according  to  specific  needs  are  more  critical  than  the 
costs. 

As    the    use    of     composite    materials     in  structural 
applications   increases,    there   is  more  need   for  structural 
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analysis.  Unlike  an  isotropic  beam,  stretching  a  laminated 
composite  beam  can  produce  bending,  shear  and  twisting  in 
addition  to  the  extensional  deformation.  The  coupling  among 
the  modes  of  deformation  is  the  major  obstacle  to  precisely 
modeling  the  mechanical  behavior  of  beam  like  composite 
structures . 

In  the  present  chapter,  a  new  and  efficient  method  to 
analyze  delaminations  in  laminated  composite  beams  is 
presented.  The  equations  of  equilibrium  were  derived  by  Sankar 
(1991)  from  the  classical  shear  deformable  laminated  plate 
theory  (Whitney,  1987).  The  equilibrium  equations  are  assumed 
to  be  satisfied  in  an  average  sense  over  the  width  of  the 
beam.  From  this  assumption  results  a  new  set  of  force  and 
moment  resultants,  which  amplify  the  possibilities  for 
modeling  beams  under  different  loading  conditions,  including 
torsion.  The  equilibrium  equations  are  derived  from  the 
Minimum  Potential  Energy  Principle. 

In  the  present  study  a  new  offset  beam  finite  element  is 
developed  for  modeling  the  problem.  Initially,  the  capability 
of  the  new  finite  element  in  modeling  the  problem  is  verified 
by  solving  the  problem  of  a  cantilever  beam  under  different 
end  loading  conditions.  The  results  are  compared  with  those 
from  beam  theory  solutions  found  in  Timoshenko  (1970), 
Reismann  and  Pawlik  (1980),  Whitney  (1987)  and  Sankar  (1991). 
Next,  the  strain  energy  release  rate  for  a  delaminated 
specially    orthotropic    cantilever    beam    under    torsion  is 
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calculated.  The  delamination  is  considered  to  be  at  the 
middle-plane  of  the  beam.  The  results  for  G  from  the  finite 
element  method  are  compared  with  analytical  solutions  of  the 
problem.  The  effect  of  position  of  delamination  on  the  strain 
energy  release  rate  values  is  presented.  Finally,  the  problems 
of  a  specially  orthotropic  and  a  quasi-isotropic  simply 
supported  delaminated  beam  under  a  transverse  force  is 
considered.  The  strain  energy  release  rate  values  for  various 
positions  of  the  delamination  are  computed  by  using  the  finite 
elements  developed  in  the  present  study.  The  results  will  be 
useful  in  explaining  delamination  propagation  in  composite 
beams  due  to  quasi-static  impact. 

3 . 2  Beam  Equations 

Considering  the  laminated  beam  shown  in  Fig.  3.1,  the 
Shear  Deformation  Theory  (Whitney,  1987)  for  laminated  plates 
assumes  the  displacement  field  as 

uix.y,  z)  -Uo(x,y)  +zilr^(x,  y) 

v{x,y,z)'Vo  ^/y) +zty(x,y)  (3.1) 
w(x,y,  z)  -Woi^.y) 

where  Ug,  Vq  and  Wg  are  the  displacements  in  the  middle  plane 
of  the  plane  in  the  x,  y  and  z  directions,  respectively,  and 
so,  are  tjj^  and  t^^  with  respect  to  the  rotations  about  the  X  and 
Y  axes. 

All  five  displacements  and  rotations  in  equation  3.1  can 
be  expanded  as  a  Taylor  series,  in  the  direction  of  the  width 
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Figure  3 . 1  Laminated  composite  beam 


of  the  plate.  If  just  the  first  order  terms  of  this  series  are 
retained,  the  displacements  and  rotations  can  be  expressed  in 
the  form  of 

Ug  {x,y)  -C7lx)  +yF{x) 
Vo(x,y)  +yG(x) 

w  U.y) -W{x) +ye{x)  (3.2) 
i|r^{x,y)  -<j)(x)  +ya  (x) 

y{x,y)  -W  ix)  +yP(x) 

The  displacement  field  is  now  represented  in  terms  of  ten 
functions  of  the  x  coordinate. 

u{x,y,z)  -C7(x)  +yF(x)  +2A{x)  +yza  (x) 
v{x,y,z)-V{x)  +yG{x)  +zT(x)  +yzP(x) 
w(x,  y,  z)  -  W{x)  +ye  (x) 

From  the  beam  theory,  we  assume  that  there  is  no 
deformation  in  cross  sections  normal  to  the  x  direction  (e^^=0, 
eyy=0  and  Kyj=0)  .  From  this  assumption  and  Eq.  3.3  results 


and 


-  Gix)  +  ZP(X)  -  0  (3.4) 

Yy;,-i|r  (x)  +yp  (x) +0  (x) -0  (3.5) 
From  equation  3.4  and  3.5  we  obtain  that 

G(x)-0 

P(x)-0  (3.6) 

Y(x)— e(x) 

The  displacement  field  can  now  be  represented  in  terms  of 
seven  functions  of  the  x  coordinate,  instead  of  the  ten 
functions  used  before. 

"(-x:,y,  z)  -f7(x)  +yF(x)  +z(j)(x)  +yza  (x) 
v{x,y,z)-V{x)  -zdix) 
w{x,y,  z)  -W{x)  +y8(x) 


Denoting  the  differentiation  with  respect  to  x  by  a 
prime,  we  represent  the  deformation  field  from  Eq.  3.7  as 


e 

XX 
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Y 

xy 

[y 

ZX) 

\ 

(j)'  +  ya' 

(^x°] 

+  z 

yxy° 

+  Z 

^  xy 

J 

[     0  j 

Y  ° 

I  0  ) 

U'^yF' 
F+  V' 


The  deformation  field  can  also  be  expressed  as 

E'-'E+yE 


where 


(e  ) 

(  U'  ] 

{f'\ 

y  xyO 

V'+F 

0 

!  (E)  - 

; 

a' 

^xy 

a-e' 

0 

^y  ZX, 

U+e'j 

(3.8) 


(3.9) 


(3.10) 


From  the  constitutive  relations  the  stresses  developed  in 
the  k**'  layer  are  given  by 


^"'\-[Q]\^''°Uz[Q] 

xy)  V  «  xyoj 


(3.11) 


and 


^zx-  Qss  (<l>+^')  +  y  c>55  (a+e') 


(3.12) 


where 


[Q]  - 


%1    Ql2  ^16^ 
^12    Q22  Q26 
Pl6    Q2S  OsS, 


(3.13) 


62 

The  terms  Q..  in  Eq.  3.13  are  stiffness  terms  and  can  be  found 
in  Whitney  (1987) . 

By  definition,  the  force  (N)  and  moment  (M)  acting  on  the 
laminate  are 


'  2 
2 


{N^,M^)- f  X^[l,z]dz  (3.14) 


2 
Ji 
2 

'  2 

and  they  are  shown  in  Fig.  3.2. 

The  stiffness  coefficient  terms  are  defined  as 


{A^j,B^j,D^j)-  f  Q,j[l,z,z^]  dz  (3.15) 


_  J3 
2 

By  using  equations  3.9-13  and  the  definition  of  stiffness 
coefficients,  the  Eq.  3.14  can  be  written  in  a  more  simplified 
form  as 

F-C{E+y^)  (3.16) 

where 

(F)  ^-  {N^,N^,M^,M^,V^)  (3.17) 
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Figure  3.2  Force  and  moment  resultants 


and  the  matrix  C  is  formed  with  the  following  stiffness 
coefficient  terms 


^16 

fill 

file 

0  ^ 

A,s 

^ee 

file 

fiee 

0 

fill 

file 

^11 

^le 

0 

file 

fiee 

^le 

^ee 

0 

1° 

0 

0 

0 

The  factor  k  presented  in  term  C55  is  the  shear  correction 
factor  (Whitney,  1987),  and  it  is  considered  equal  to  5/6  in 
the  present  work. 


3.3  New  Set  of  Force  and  Moment  Resultants 

A  new  set  of  force  and  moment  resultants  is  defined 
from  the  integration  of  the  column  vector  of  forces  along  the 
width  of  the  beam: 

2 

F-  I  Fdy  (3.19) 

2 


and 


b 

2 


(3.20) 


In  the  two  last  equations,  the  vector  of  forces  F  is 
given  by  equation  3.17,  and  b  is  the  width  of  the  beam. 

Substituting  equation  3,16  in  Eq.  3.19-3.20,  and  assuming 
that  the  stiffness  C  of  the  beam  is  constant,  we  obtain 


65 


F'  f  C  Edy-bCE 


-b/2 
biz 

F-  f  C  E  y^dy 
-i/2 


•(«) 


CE 


The  explicit  beam  contitutive  relations  are: 


(3.21) 


N 


[F]- 


xy 


M 


xy 


V 
\  xi 


-b 


■^11   ■^16    -^11  -^16 


■^16    -^6    ^16  ^66 

•^ii  -Bie  -Oil  ^16 
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0  icM, 
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a-Q' 


(3.22) 
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\  X/ 
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■^11  -^16  ^11  ^16  0 

■^16  -^6  ^16  ^66  0 

^11  ^16  ^11  ^16  0 

^16  ^66  ^16  ^66  0 


0  k^A, 


55/ 


'  F' 
a' 
0 
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(3.23) 


3.4  Equilibrium  Equations 

The  Principle  of  Minimum  Potential  Energy  is  applied  in 
the  derivation  of  the  equilibrium  equations.  The  principle 
says  that  of  all  displacement  fields  which  satisfy  the 
prescribed  constraint  conditions,  the  correct  state  is  that 
which  makes  the  total  potential  energy  of  the  structure  a 
minimum  (Tauchert,  1970) .  The  total  potential  energy  of  the 
structure  (n)  is  obtained  from  the  sum  of  the  strain  energy 
of  the  beam  (0)  and  the  potential  of  the  external  force  (/)  : 
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71  -  (|)  +  X 


(3.24) 


According  to  the  Principle  of  Minimum  Potential  Energy, 


where  the  symbol  S  is  the  variational  operator. 

Considering  only  transverse  loading  acting  on  the  beam 
surface  in  the  z-direction,  the  potential  relative  to  the 
external  loading  is  given  by 


where  q(x,y)  and  w(x,y)  are  the  transverse  loading  and 
displacement,  respectively,  and  L  is  the  beam  length. 

Similar  to  the  procedure  adopted  to  define  a  new  set  of 
force  resultants  in  Eqs  3.19-20,  the  transverse  loading  is 
divided  into  two  parts  and  defined  by 


Sti  -  6<|>  +  6x  -  0 


(3.25) 


b 

L  2 


(3.26) 


0 

2 


Qix,y)  -  g{x)  +  y  <S(x) 


(3.27) 


where 


2 


2 
b 

2 


(3.28) 
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Substituting  the  definition  of  transverse  displacement 
(Eq.  3.7)  and  transverse  loading  (Eq.  3.28)  in  Eq.  3.26  gives 
the  result 


X  -  -f{g(x)W(x)  +  (S(x)QU))  dx  (3.29) 


From  Eq.   3.16,   the  strain  energy  per  unit  area   (0^)  of 
the  beam  can  be  derived  as 

<1)^  -  ^  E^CE 

T    _  _  (3.30) 

-  ^  (E+yE)  '^CiE+yE) 

The  strain  energy  per  unit  length   (0^)    is  obtained  by 
integrating  <p^  along  the  width  (b)  of  the  beam 

i 
2 


♦  l  -  f^A  dy 


2 

b 
2 


(3.31) 


-  -|  /  (E+yi)  ^C{E+y^  dy 


2 

2 


Because 


2 


j  E'^CEydy-0 

b  (3.32) 

2 

j  E'^CEydy-0 
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the  Eq.  3.31  takes  the  form  of 


2 


2 


4( 


2 


(3.33) 


bE'^CE  +  E^C^I 


The  total  strain  energy  in  the  beam  is  expressed  as 


<|>  -  j^L^ 


2Jl  12  ) 


(3.34) 


dx 


The  total  potential  of  energy  in  the  whole  structure  (tt) 
is  obtained  by  substituting  the  Eq.  3.29  and  3.34  into 
equation  3.24. 


bE'^CE+  —  e'^CE-  q{x)  W{x)  -  ^(x)  Q  ix) 
X  2 


dx  (3.35) 


Applying  the  Principle  of  Minimum  Potential  Energy  we 
obtain  the  following  equilibrivim  equations 
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(3.36) 
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Next,  from  each  equilibrium  equation  in  Eq.  3.36  results 
one  corresponding  natural  boundary  condition,  which  is  shown 
in  the  next  equations  as 


dx 

-0 

N^bU-0 

bx 

-0 

N^bV-0 

dx 

—  g 

V^bW-0 

dx 

N^bF-0 

dx 

(V,-M^)  66-0 

dM, 

dx 

M^6(j)-0 

dM, 
dx 

M^ba-0 

(3.37) 


By  using  Eqs.  3.22-23,  the  new  set  of  force  and  moment 
resultants  can  be  expressed  in  terms  of  seven  unknown 
functions:  U,  V,  W,  F,  <p,  a  and  6.  Substituting  for  the 
resultants  in  terms  of  displacements  in  the  equilibrium 
equations  3.37  results  in  a  system  of  seven  ordinary 
differential  equations  for  the  displacements. 

3 . 5  Finite  Element  Solution 

The  system  of  ordinary  differential  equations  obtained 
from  the  equilibrium  equations  3.37  is  solved  in  this  work 
using  the  finite  element  method.  A  new  finite  element  has 
been  developed  to  model  the  beam  problem.  This  finite  element 
has  three  nodes  and  seven  degrees  of  freedom  (U,  V,  W,  F,  <p, 
a,  e)   at  each  node.  The  middle  node  is  statically  condensed 
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when  solving  the  problem  for  displacements,  but  it  is 
considered  when  calculating  the  strain  energy  release  rate,  in 
order  to  obtain  a  more  accurate  solution.  The  nodal  forces  and 
moments  for  the  i*''  node  of  the  structure  are  F^^,  F^,  F^,  M^,  M^, 
W  and  T. 

A   quadratic    variation    of    all    seven    displacements  is 
assumed  along  the  element  length.  Denoting  by  X  any  specific 
displacement,  and  by  X.  this  specific  displacement  at  the  i*'' 
node,  the  displacement  and  deformation  within  the  element  are 
defined  as 

x'-b^Xj^-^b^X^^b^X^  IJ.je; 
In    the    last    equation,    the    terms    a,,  and    aj  are 

interpolation  functions  in  the  form  of 

a^-l-3x+2}? 

a^-x{2x-l)  (3.39) 
a^-4.x{l-2x) 

and  b^,  and  b3  are  their  derivatives  with  respect  to  the  x 
direction 


i^i-  (4x-3)/L 
jbj-  (4X-1)/L 
2)3-4  (l-2x)/L 


(3.40) 


71 


where 

In  the  last  equation,  L  is  the  element  length  and  x  is 
the  coordinate  along  the  beam  axis. 

By  using  Eq.  3.38,  and  the  terms  a,  and  b,.  (i=l,2,3) 
defined  in  Eqs  3.39  and  3.40,  the  deformation  field  (Eq.  3.10) 
within  each  finite  element  as  a  function  of  the  nodal 
displacements  is  expressed  in  the  form  of 

(e)  -  [a]  (g)  ,3  v 

where  q  is  the  vector  of  the  displacements 

U^.V^.W^,F^,^^,a^,Q^,     ,  (3.42) 

and  a  and  a  are  two  strain-displacement  matrices,  both  with 
dimension  of  (5  x  21)  . 

The  first  strain-displacement  matrix,  a,  is  composed  by 
the  following  terms: 

00000     0  252   00000  0^)3   00000  0^ 

0  2?,   Oa^OO     0  0         0   32   00  0  Oi^jOajOO  0 

0000   23iO     0  0000   2)2   0  0  0000   233    0  0 

0    0    0    0    0   a,  -2?,  0    0    0    0    0    a2  -2?2  0    0    0    0    0   53  -2)3 

00   2?i0ai   0     0  00   2^2    0    32    0  0  00   2^3    0   33    0  0 

(3.43) 
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and  the  second  one,  a,  of: 

'OOOi?iOOO  OOOijjOOO  OOOJbjOOON 

000000    0  000000    0  000000  0 

00000  jb^O  00000  i^gO  00000  jb3  0 

000000    0  000000    0  000000  0 

0  0  0  0  0       i?^  0  0  0  0  0  a2  0  0  0  0  0  2^3^ 

(3.44) 

similarly  to  Eq.  3.10,  the  total  deformation  within  each 
element  is  given  by  the  expression 


e-e+yS 


(3.45) 


or,  after  using  the  strain-displacement  matrices,  by 


e-  (a  +  yi)  q 


(3.46) 


With 


a-a+yi 


(3.47) 


the  element  stiffness  matrix  is  defined  as 
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In  Eq.  3.48,  D  is  the  elasticity  matrix  for  composite 
materials.  The  terms  in  matrix  D  are  the  same  terms  shown  in 
Eq.   3.18.  Each  term  of  the  matrix  D  is  defined  as 

b 

{Aij,B^j,D^j)  -Jq^j  [1,Z,z2]  dz  (3.49) 
a 

where  a  and  b  are  respectively  the  lowest  and  hiqhest 
coordinates  in  the  z  direction  in  the  considered  element.  By 
defining  two  element  stiffness  matrices     in  the  form  of 

L 

K'bfa'^D  a  dx 

°  ^  (3.50) 

K-^  (S-'D  a  dx 
12  J 

0 

the  equation  3.48  is  reduced  to 

K  -K  +  K  (3.51) 

The  elements  of  the  stiffness  matrices  K  and  K  are  shown 
in  the  Appendix.  The  dimension  of  the  stiffness  matrix  k  are 
reduced  from  (21  x  21)  to  (14  x  14)  by  using  static 
condensation  and  the  stiffness  matrix  of  each  finite  element 
is  assembled  in  a  global  stiffness  matrix,  representative  of 
the  stiffness  of  the  whole  structure. 

The  final  representation  of  the  problem  is 

F  -  K  q  (3.52) 


74 


where  K  is  the  total  stiffness  matrix,  and  F  and  q  are  vectors 
with  the  nodal  forces  and  displacements,  respectively.  In  the 
present  work,  the  equation  3.52  is  solved  using  the  Gauss 
Elimination  Method  (Bathe,  1982) . 

3.6  Strain  Energy  Release  Rate 

The  strain  energy  release  rate  is  calculated  by  using  the 
expression  for  J-integral.  A  zero  volume  path  is  delineated 
surrounding  the  crack  tip,  as  shown  in  Fig.  2.8,  and  the  J- 
integral  expressed  in  terms  of  the  strain  energy  per  unit 
length  (U)  of  each  element  around  the  crack,  as  in  Eq.  2.55, 


The  strain  energy  per  unit  length  of  the  beam  (Eq.  3.34) 
is  given  by 


The  expression  for  the  strain  energy  per  unit  length  of 
the  i*"  element  is  obtained  by  substituting  Eq.  3.46  in  the 
previous  equation. 


(3.53) 


(3.54) 


(3.55) 


In  Eq.  3.55,  the  vector  q  contains  the  displacements  of 
the  nodes  of  the  i^*^  element,  including  the  middle  node.  The 
reason  for  considering  the  middle  node  in  the  calculation  of 
the  strain  energy  is  because  the  static  condensation  process 
assumes  that  no  forces  or  couples  are  applied  to  the 
intermediate  node.  This  assumption  is  not  true  for  the  case 
where  bottom  and  top  offset  elements  are  put  together  to  model 
the  beam. 

3 . 7  Rigid  and  Gap  Elements 

The  possibility  of  interference  between  the  crack 
surfaces  in  the  delaminated  region,  when  the  beam  is  under 
loading,  principally  under  torsion  loading,  requires  that  gap 
elements  or  rigid  elements  be  placed  in  appropriate  positions 
in  order  to  monitor  the  contact. 

Two  types  of  rigid  elements  are  anticipated  to  be  used 
depending  on  which  side,  left  or  right,  of  the  beam  presents 
the  interference  phenomenon.  Figure  3.3  illustrates  the 
interference  situation  and  the  rigid  element  actions. 

In  matrix  form,  the  equations  of  the  rigid  element  placed 
on  the  left  side  can  be  represented  as 
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where  t  and  b  represent  top  and  bottom  elements,  respectively, 
and  the  force  F  is  shown  in  Fig.  3.3. 

Similarly,  the  equations  of  the  rigid  element  placed  to 
on  the  right  side  are 
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The  gap  element  is  placed  in  the  structure  to  avoid  the 
interference  in  the  delaminated  region  while  the  structure  is 
under  vertical  loading.  In  matrix  form,  the  constraint 
equations  relative  to  the  gap  element  are: 
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1  ^ 
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(F.)b 
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-1 
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(3.58) 
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Figure  3.3  Rigid  elements:   (a)  Beam  under  torsion 

(b)  Interference  on  left  side  (front  view) 

(c)  Interference  on  right  side  (front  view) 
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The  rigid  element  and  gap  element  matrices  in  Eqs.  3.56- 
3.58  are  assembled  in  the  global  stiffness  matrix  using  the 
same  procedure  adopted  to  assemble  the  element  stiffness 
matrices.  In  Eqs.  3.56-3.58  is  the  contact  force  between 
the  interfering  elements,  or  is  the  couple  due  to  the 
contact  force  and  F  is  the  force  in  the  rigid  or  gap  element 
in  the  z  direction. 


CHAPTER  4 
NUMERICAL  RESULTS  AND  DISCUSSION 


4.1  Numerical  Results 

This  chapter  presents  the  numerical  results  obtained  from 
the  finite  element  described  in  the  previous  chapter.  Three 
different  tests  were  performed,  and  their  numerical  results 
compared  with  those  obtained  from  the  beam  theory. 

The  first  test  of  the  beam  finite  element  consists  of  a 
specially  orthotropic  cantilever  beam  subjected  to  several 
different  end  loading  conditions.  This  test  is  divided  into 
two  parts:  one  in  which  the  beam  is  modeled  using  two 
sublaminates  (top  and  bottom)  and  the  other,  where  the  beam  is 
formed  by  just  the  top  sublaminate.  The  purpose  of  theses 
tests  is  to  prove  the  validity  of  the  finite  element  developed 
in  modeling  the  beam  problem,  including  the  torsion  problem. 

In  the  second  numerical  test,  a  specially  orthotropic 
delaminated  cantilever  beam  is  put  under  torsion  loading. 
Initially,  the  delamination  is  supposed  to  be  in  the  middle 
plane  of  the  beam,  and  the  result  in  terms  of  the  strain 
energy  release  rate  is  compared  with  the  closed-form  solution. 
The  delamination  is  then  placed  between  different  layers  along 
the  thickness,    and  a   study   is  performed  to  understand  the 
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effect  of  the  delamination  position  on  the  strain  energy 
release  rate. 

The  last  numerical  test  consists  in  the  simulation  of  a 
quasi-static  impact  test  in  specially  orthotropic  and  quasi- 
isotropic  simply  supported  beams.  Some  experimental 
observations  for  this  specific  tests  are  available  in  our 
Composite  Laboratory  (Kwon,  1991) ,  and  the  present  numerical 
study  will  help  interpret  the  experimental  results.  Again,  the 
focus  was  on  the  effect  of  the  delamination  position  on  the 
values  of  the  strain  energy  release  rate. 

4.2  Material  Properties  and  Stiffness  Coefficients 

In  all  the  numerical  tests  performed  in  the  present 
study,  the  beam  is  supposed  to  be  made  up  of  material  with  the 
following  properties: 

Longitudinal  Elastic  Modulus  (E^)  :  14.00  GP^ 

Transverse  Elastic  Modulus       (Ej)  :  1.00  GP^ 

Shear  Modulus  (Gij)  :  0.53  GP^ 

Poisson's  Ratio  (v^^) o.30 

Poisson's  Ratio  (v^^) o.55 

The  above  values  are  typical  values  for  high  performance 
graphite/ epoxy  unidirectional  composites. 

Table  4.1  lists  the  laminate  stiffness  coefficients  for 
a  specially  orthotropic  beam  with  thickness  equal  to  2  mm, 
and  the  reference  X-Y  plane  considered  to  be  in  the  middle 


81 

(MIDDLE)  ,  on  the  top  (BOTTOM)  or  on  the  bottom  (TOP)  surfaces 
of  the  beam. 

4.3  Cantilever  Beam 

The  dimensions  of  the  beam  are  0.2  m  x  0.02  m  x  0.002  m 
(length  x  width  x  thickness) ,  and  the  beam  is  subjected  to  a 
set  of  six  different  end  loading  conditions,  as  shown  in  Fig. 
4.1  and  Fig.  4.3. 

4.3.1  Model  with  two  sublaminates 

Initially,  the  beam  subjected  to  the  forces  shown  in  Fig. 
4.1  is  modeled  by  using  two  connected  sublaminates  (top  and 
bottom),  as  shown  in  Fig.  4.2.  Table  4.2  presents  the 
numerical  and  beam  theory  solutions,  in  terms  of  the 
displacements  at  the  tip  of  the  beam  for  each  of  the  loading 
conditions. 

The  expressions  for  beam  theory  results  can  be  found  in 
Timoshenko  (1970),  Reissman  and  Pawlik  (1980)  and  Sankar 
(1991)  ,  or  can  be  derived  from  the  beam  equations  presented  in 
Sankar  (1991) .  In  the  following  expressions  the  superscripts 
t  and  b  stands  for  top  and  bottom,  respectively. 

4.3.1.1  Beam  under  Force  1  fF^) 

In  this  case,  we  have  just  the  U  displacement.  All  other 
displacements  are  zero. 
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Table  4.1  Stiffness  coefficients 


1  STIFFNESS 
COEFFICIENTS 

POSITION  OF  THE  X  AXIS 

TOP 

BOTTOM 

MIDDLE 

Ai1 

0.2818173E8 

0.2818173E8 

0.2818173E8 

^6 

0.0 

0.0 

0.0 

0. 1066621E7 

0. 1066621E7 

0. 1066621E7 

Bii 

0.2818229E5 

-0.2818229E5 

0.0 

16 

0.0 

0.0 

0.0 

0. 1066643E4 

-0. 1066643E4 

0.0 

D„ 

0.3757714E2 

0.3757714E2 

0,9394285E1 

0.0 

0.0 

0.0 

0. 1422219E1 

0. 1422219E1 

0.3555547 

K2  A55 

0. 8888511E6 

0.8888511E6 

0.8888511E6 
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Figure  4.1  Loading  conditions  at  middle  plane 
(a) Perspective  view     (b) Front  view 
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Figure  4.2  Model  with  42  elements 
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L  (4.1) 


b  {A  \^^A 


4.3.1,2  Beam  under  Force  2  (F..) 

Under  Force  2  the  beam  has  the  displacements  V  and  F.  The 
remaining  five  displacements  are  zero. 


bHA\^^A\^)  b{A'=,^^A^,^) 

(4.3) 


F--F, 


b^{A  \^*A  \^) 


4.3.1.3  Beam  under  Force  3  fF,) 

From  this  case  results  the  displacements  W  and  <p. 

W-fJ  .  2L  \  (4.4) 


3b(D\,^D\,)  b{A%,^A\,)  j 
2  2j(i?t^^+z?*^^) 


4.3.1.4  Beam  under  Force  4  (F^) 

The    displacements    under    Force    4    are    V    and    F.  The 
remaining  displacements  are  zero. 
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12L 


(4.7) 


4.3.1.5  Beam  under  Force  5  fF^) 

Under  this  loading  condition  the  displacements  are  W  and 
0,  and  all  others  are  zero. 


1  7"  2 


(4.8) 


2  b(D\,*D\,) 


(4.9) 


<t)-F5 


4.3.1.6  Beam  under  Force  6  (F^)  -  " 

The    displacements    under    Force    6    are    6    and    a.  The 
displacement  dis  given  by  the  expression  (Sankar,  1991b) 


4b{D%,^D\,) 


(4.10) 


where 


(4.11) 


and 


while  the  displacement  a  is  defined  as 


12 


o-C,cosh(x)  +C,sinh(x)  -C,  


where 


12   55 


C2 — Cj^A,sinh(;<) 


and 


In  Eqs.  4.13-4.16, 
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(4.19) 


(4.20) 


D 


•* 


66 


4.3.2  Model  with  top  sublaminate 

In  a  second  test,  the  beam  is  modeled  using  top  offset 
elements  only  (Fig.  4.3).  The  beam  is  subjected  to  the  set  of 
loading  conditions  displayed  in  Fig.  4.3.  Table  4.3  shows  the 
numerical  and  closed  solutions  for  this  test.  The  sources  for 
the  beam  theory  solutions  presented  in  Table  4 . 3  are  the  same 
used  before:  Timoshenko  (1970),  Reismann  and  Pawlik  (1980)  and 
Sankar  (1991b) .  Some  expressions  are  derived  from  the  beam 
equations. 

4.3.2.1  Beam  under  Force  1  (F^) 

The    tip  displacements  under  F,  are  U,  W  and  0.  The 
remaining  displacements  are  zero. 


L 


(4.21) 


11 


W  -  ^ 


(4.22) 


Figure  4.3  Loading  conditions  at  the  bottom  plane 
(a)  Perspective  view  (b)  Front  view 
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cantilever  beam 


Z  A 


top  offset  elements 


Figure  4.4  Model  with  21  elements 
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(|)  L 


(4.23) 


4.3.2.2  Beam  under  Force  2  (F.) 

In  the  case  of  a  cantilever  beam  under  Fj  the  tip 
displacements  are  V,  F,  ^  and  a. 

The  displacement  V  can  be  represented  as 


V-  V^^ 


(4.24) 


where 


V,  -  (F2) 


4L3     ^     2L  \ 


3  a  c 


n 


55  J 


(4.25) 


and 


(4.26) 


The  displacement  0  is  necessary  to  calculate        and  is 


defined  as 


e  -  F, 


2 
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1  + 


55  / 


where  the  term  D^jj  is  defined  as 


(4.27) 


55 


12 


55 


(4.28) 


The  displacement  F  is  given  by 
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6L^ 
3  a  t 


(4.29) 


11 


Finally,  the  expression  for  the  displacement  a  is 


a  -  qcosh(x)  +C2sinh(z)  -C3 


48 


(4.30) 


where 


C'l-n 


48  D' 


55 


(4.31) 


C2  -  -qA,sinh(^) 


(4.32) 


and 


55 


66 


In  Eqs.  4.30-4.33, 


(4.33) 


X  '  L 


48  D\eD\^ 


(4.34) 


12£>' 


66 


11  / 


(4.35) 


55 


12 


66 


(4.36) 
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D 


•  * 


66 


66 


(4.37) 


and  the  term  D%,  is  defined  as 


55 


(4.38) 


4.3.2.3  Beam  under  Force  3  (F,) 


Under  force  Fj,  the  beam  presents  three  displacements:  U, 
W  and  0.  They  are  defined  as 


4.3.2.4  Beam  under  Force  4  CF^) 

The  displacements  in  this  case  are  V,  F,  a  and  6.  The 
remaining  displacements  are  zero.  This  case  can  be  considered 
as  the  result  of  the  addition  of  two  different  loadings,  I  and 
II,  as  shown  in  Figure  4.5. 


(4.39) 


and 


2 


(4.41) 
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h/2b 


The  displacements  for  the  Loading  I  were  calculated  as 
explained  in  the  Section  4.3.1.4.  From  Loading  II  results  the 
displacements  a  and  8,  which  are  calculated  using  the 
following  expressions  (Sankar,  1991b) : 

(4.42) 

a-C8sinh(A.L) 

where  Cg  is  a  constant  defined  as 

^      12/?„  1  (4.43) 

^    jb^n  Xcosh(AL) 


1 
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where  the  constant  A  is  calculated  using  Eq.  4.18. 

From  the  definition  given  in  Eq.  3.20,  we  have  that 


2 

M,-  f  M,dy 

-b/2 

-2^4 
2  2b 

h 


(4.44) 


The  rotation  6  is  given  by 


where  the  constant  Cj  is  defined  in  Eq.  4.16. 

Considering  the  boundary  condition  e(0)=0,  we  have 

e(x)  —c^aL 

Loading  I  does  not  produce  a  or  e  displacements.  Hence, 
the  a  and  9  displacements  calculated  using  the  previous 
equations  are  the  result  for  the  total  loading. 

The  displacement  F  is  given  by  Eq.  4.7  and  the 
displacement  V  is  obtained  using  the  expression 

v-v,  +6  — 
where  V,  is  given  by  Eq.  4.6. 

4.3.2.5  Beam  under  Force  5  (F-) 

Under  force  Fg  the  displacements  are  U,    W  and  0.  The 
remaining  displacements  are  zero, 
where  the  displacement  <p  is  given  by 
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U  -  -^<1)  (4.48) 


♦-fs^^  (4.49) 


The  displacement  W  is  defined  as 


^--^.-H^  (4.50) 


4.3.2.6  Beam  under  Force  6  (F^) 

In  this  case  the  displacements  are  V,  F,  a  and  d. 


2 


(4.51) 


Where  the  rotation  6  is  calculated  using  the  expression 


e  - 


b  D 


66 


1  + 


66 


4  D 


55  / 


(4.52) 


The  rotation  F  is  defined  as 


F  -  tan-1  (  — ) 
L 


(4.53) 


and  the  warping  function  a  can  be  obtained  using  Eqs.  4.30- 
4.37. 


4.3.3  Comments  on  the  finite  element  so Tut inn 

Tables  4.2  and  4.3  show  a  very  good  agreement  between  the 
present  finite  element  solutions  and  the  analytical  results 
derived  from  beam  theory.  This  agreement  can  also  be  seen  in 
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Figures  4.6  and  4.7,  in  which  is  shown  the  angle  of  twist 
along  the  beam  length  for  two  cantilever  beams  with  different 
lengths.  The  beams  are  under  a  unit  torque,  and  in  both  cases 
they  were  modeled  by  ten  finite  elements.  From  these  figures, 
we  can  see  that  in  the  region  close  to  the  support  the 
agreement  is  not  good,  and  probably  it  will  require  use  of 
larger  number  of  finite  elements  near  the  fixed  support. 

The  beam  theory  solution  for  this  case  (Sankar,  1991b) 
incorporates  the  edge  effect,  represented  by  the  third  term  in 
brackets  on  the  right  hand  side  of  Eq.  4.52,  and  is  given  by 
the  following  expression 


where  e(x)  and  e(L)  are  the  angle  of  twist  at  a  distance  x 
from  the  fixed  support  and  at  the  tip  of  the  beam, 
respectively. 

The  tip  rotation  of  a  beam,  e(L),  is  given  by 


e  ix)  -F, 


X 


QiL) 


(4.54) 


4  2?  A 


66 


(4.55) 


In  Eqs.  4.54  and  4.55  we  have  that 


(4.56) 


X  - 


100 


101 


M  *  «  N  T-  O 

O  O  O  O  O 

o  o  d  d  o 

(pej)e  'UOUB^OU 
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and 


55 


12 


55 


(4.57) 


^66   -  ^66  + 


12 


55 


where  is  a  shear  correction  factor,  considered  equal  to  5/6 
in  the  present  work. 


4.4  Delaminated  Cantilever  Beam  under  Torsion 

In  the  present  application,  the  beam  is  assumed  to  have 
a  delamination  of  length  20  mm.  The  dimensions  of  the  beam 
are  0.60  m  x  0.02  m  x  0.002  m.  The  length  of  the  beam  and  of 
the  delamination  were  enlarged  if  compared  with  the  beam  used 
in  the  previous  examples.  The  reason  was  to  avoid  the  edge 
effect,  so  as  to  be  able  to  model  the  problem  using  fewer 
finite  elements. 

A  unit  torque  is  applied  at  the  free  end  of  the  beam,  and 
the  strain  energy  released  is  computed.  The  numerical  result 
is  then  compared  with  the  closed-form  solution  obtained  using 
the  beam  theory  (Sankar,  1991b).  Figure  4.8,  (a)  and  (b)  , 
shows  the  beam  under  torsion  and  the  model  used  to  obtain  the 
closed-form  solution,  while  Fig.  4.9  illustrates  the  finite 
element  model  used  to  solve  the  problem.  In  this  model, 
middle-plane  beam  finite  elements  are  used  in  the  uncracked 
regions,  behind  and  beyond  the  delamination,  (regions  1  and 
3)  ,  while  top  and  bottom  offset  elements  are  used  in  the 
cracked     portion     of     the     beam     (region     2).     Later,  the 


103 


Figure  4.8  Delaminated  beam 

(a)  Under  torsion     (b) Separated  by  parts 
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Figure  4.9  Finite  element  model 
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delamination  will  be  assumed  in  different  positions  along  the 
thickness  of  the  beam,  the  strain  energy  release  rate  will  be 
computed  for  each  of  these  positions. 

4.4.1  Delamination  at  the  middle-plane  of  the  beam 

The  finite  element  model  used  to  solve  this  case  has 
forty-six  nodes  and  forty-six  elements,  as  shown  in  Fig.  4.9. 
There  are  twenty-two  regular  finite  elements  and  twenty-four 
offset  elements  (twelve  top  elements  and  twelve  bottom 
elements) . 

Referring   to   Fig.    4.8(b)    and   to   the      strain  energy 
release  rate  (G)  definition  (  Hellan,  1984)  we  have 

~  dA~  b  dL^ 

where  U  is  the  strain  energy  of  the  beam  and  A  is  the  crack 
area.  In  the  present  case,  the  area  of  the  crack  is  formed  by 
a  rectangle  with  dimensions  egual  to  b  and  Lj. 

The  strain  energy  of  the  beam  (U)   is  defined  as 

U-  ^Te(L)  (4.59) 

where  T  is  the  torque  and  5(L)  is  the  rotation  e  at  the  tip  of 
the  beam. 

The  rotation  at  the  tip  of  the  beam  is  calculated  by 
adding  the  rotations  in  the  three  parts  of  the  beam  shown  in 
Figure  4.7,  and  is  defined  by 
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(4.60) 


Where  the  subscripts  stand  for  the  regions  in  the  beam  (Fig. 
4.8(b) ) . 

By  definition, 


the  superior  and  inferior  limits  of  the  beam  in  the  z 
direction. 

The  parts  l  and  3  of  the  beam  have  the  same  integration 
limits,  h/2  and  -h/2,  while  for  part  2  they  are  h/4  and  -h/4. 
Considering  these  integration  limits  and  solving  Eq.  4.61,  we 
obtain  the  following  relations  among  the  stiffness 
coefficients  of  the  three  parts  of  the  beam: 


By  using  Eqs.   4.54,   4.55  and  4.60,   and  neglecting  the 
edge   effect,    the   rotation   at   the   tip   of   the   beam   is  now 


(4.61) 


where  the  integration  limits  h^  and  h.  stand  respectively  for 


(4.62) 


expressed  as 


e(L) 


 n    +   + 


4b      [d,,  D, 


1 


(4.63) 


4  A 


55 


Referring  to  Fig.  4.8(b),  we  have  that 
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L1+L3-L-L2 

Substituting  Eq.  4.64  in  Eq.  4.63  gives 

r(L-L,)  I 


(4. 64) 


e(L)  - 


4i) 


1  .  1 


(4.65) 


Using  Eq.4.59  we  obtain  the  strain  energy     of  the  beam 
(U)  as 


+  

{  ^ 

1^66        ^55  j 

b 

8Z?. 


(4.66) 


55/ 


Considering  Eq.  4.58,  we  obtain  the  following  expression 
for  strain  energy  release  rate  (G) 


G  - 


rp2 


1 


D, 


55/ 


2^6  ^D^^i, 


(4.67) 


or,  after  simplifying  the  last  equation, 

3  r2 


G  - 


8  b^n 


(4.68) 


66 


Substituting  in  the  last  equation  the  values  of 


r-1  Nm 
b-0. 02m 


Z?66  -0-3556  -^m 


N  _3 
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we  obtain  for  the  strain  energy  release  rate  the  value  of 
2,636.0  J/m^,  while  the  numerical  result  for  the  present  case 
is  2,545.5  J/m^.  The  relative  difference  between  numerical  and 
closed-form  solutions  is  3.45%. 

4.4.2  Delamination  position  varying  along  the  thickness 

In  this  example,  the  delamination  was  assumed  to  be  in 
different  positions  along  the  thickness  of  the  beam.  The 
procedure  adopted  here  is  similar  to  that  used  to  calculate 
the  strain  energy  release  rate  with  the  delamination  placed  in 
the  middle  plane,  except  for  the  integration  limits  when 
calculating  the  stiffness  coefficients  A..,  B.^.  and  D... 

Figure  4.10  illustrates  the  adopted  model  for  the  case  of 
delamination  placed  between  layers  close  to  the  top.  Table  4.4 
presents  the  obtained  results  for  different  positions  of  the 
delamination  along  the  thickness  of  the  beam.  These  results 
are  also  depicted  in  Figure  4.11.  From  Figure  4.11,  we  can 
conclude  that  the  possibility  of  crack  propagation  will  be 
increased  as  the  crack  position  approaches  the  middle  plane  of 
the  beam.  As  the  crack  plane  approaches  the  top  or  bottom 
faces  of  the  beam,  the  strain  energy  release  rate  becomes 
negative.  This  could  be  due  to  errors  in  the  energy 
calculations,  and  can  be  avoided  by  having  shorter  elements 
near  the  crack  tip.  This  problem  should  be  subjected  to 
further  investigations . 
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B  :  bottom  offset  element 


Figure  4.10  Delamination  close  to  the  top  plane 
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Table  4.4  Strain  energy  release  rate  (G)  for  cantilever 
beam 


POSITION 

G  (J /m  } 

between  layers  1  and  2 

-17.5 

between  layers  2  and  3 

173.2 

between  layers  3  and  4 

653.8 

between  layers  4  and  5 

1,409.9 

between  layers  5  and  6 

2,195.1 

between  layers  6  and  7 

2,545.5 

between  layers  7  and  8 

2,195.5 

between  layers  8  and  9 

1,409.9 

between  layers  9  and  10 

653.8 

between  layers  10  and  11 

173.2 

between  layers  11  and  12 

-17.5 
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4.5  Simply  Supported  Beam  under  Transverse  Force 

In  this  section  we  consider  the  problem  of  transverse 
loading  on  a  simply  supported  beam  (Figure  4.12). 

The  dimensions  of  the  beam,  the  length  of  the 
delamination,  the  material  and  the  finite  element  model  are 
the  same  as  in  the  preceding  example.  A  force  equal  to  100  N 
is  applied  to  the  center  of  the  beam.  The  finite  element  model 
is  shown  in  Figure  4.12. 

Two  classes  of  laminates  were  analyzed,  a  specially 
orthotropic  laminate  and  a  [0/45/-45/90]s  laminate.  The 
delamination  was  assumed  in  different  positions  along  the 
laminate  thickness,  and  the  results  for  each  class  of  laminate 
are  listed  in  Table  4.5  and  shown  in  Figure  4.13. 

Similarly  to  the  cantilever  beam  case,  from  the  present 
results  we  can  conclude  that  the  possibility  of  crack 
propagation  will  be  increased  as  the  crack  position  approaches 
the  middle-plane  of  the  beam.  This  phenomenon  is  less 
pronounced  in  the  specially  orthotropic  beam  than  in  the 
other  beam.  Another  conclusion  is  that  for  a  given  transverse 
force,  the  possibility  of  crack  propagation  is  more  in  the 
[0/45/-45/90]s  laminate. 

These  results  have  an  important  application  in 
understanding  the  propagation  of  delamination  damage  in 
composite  laminates  due  to  low-velocity  impact  as  well  as 
quasi-static  indentation  types  of  loading  (Kwon,   1991) . 
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R  :  regular  element 
T  :  top  offset  element 
B  :  bottom  offset  element 


Figure  4.12  Simply  supported  beam 
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Table  4.5  Strain  energy  release  rate  (G)   for  simply 
supported  beam 


1  POSITION 

G  (J/m^) 

Specially 
orthotropic 

[0/45/-45/90]s 

between  layers  1  and  2 

406.9 

1,945.0 

between  layers  2  and  3 

1,071.0 

5,858.0 

between  layers  3  and  4 

1,971.0 

10,070. 0 

between  layers  4  and  5 

2,495.0 

11,390.0 

between  layers  5  and  6 

1,971.0 

10,070.0 

between  layers  6  and  7 

1,071.0 

5,858.0 

between  layers  7  and  8 

406.9 

1,945.0 
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CHAPTER  5 
SUMMARY  AND  FUTURE  RESEARCH 

5 » 1  Summary 

The  present  investigation  was  restricted  to  delaminations 
in  laminated  beams.  Finite  elements  were  developed  to  model 
specific  cases  of  delaminations,  and  the  numerical  results 
were  compared  with  results  from  other  methods  and  approaches. 

A  laminated  shear  deformable  beam  finite  element  with 
nodes  offset  to  either  top  or  bottom  has  been  developed  to 
analyze  free  edge  delaminations.  These  offset  elements  were 
used  to  model  the  sublaminates  above  and  below  the 
delamination  plane.  A  rigid  element  was  used  to  connect  the 
nodes  at  the  crack  tip,  and  the  forces  and  moments  transmitted 
by  this  rigid  element  are  part  of  the  finite  element  solution 
of  the  problem.  The  strain  energy  release  rate  was  calculated 
using  expressions  derived  from  the  J-integral,  the  crack  tip 
force  method  and  crack  closure  method.  These  expressions  are 
in  terms  of  displacements  and  forces  in  the  elements 
surrounding  the  crack-tip,  which  are  part  of  the  finite 
element  solution.  The  present  method  considers  the  effect 
caused  by  the  thermal  stresses  induced  in  the  composite  layers 
during  the  curing  process. 
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There  are  some  advantages  in  using  the  present  method. 
The  finite  element  method  has  great  flexibility  in  modeling 
any  kind  of  geometry,  loading  and  boundary  conditions.  This 
makes  it  easier  to  model  common  specimens  used  in  laboratory 
experiments  and  compare  results.  In  modeling  the  cracked 
region  of  the  beam  by  using  offset  node  elements,  one  can 
avoid  renumbering  the  nodes  as  the  crack  propagates.  It 
tremendously  saves  computer  storage  and  time  involved  in 
redefining  nodes. 

The  delamination  in  anisotropic  beams  was  analyzed  using 
a    new    beam    finite    element    to    model    the    problem.  The 
formulation  of  this  new  beam  finite  element  is  based  on  a  beam 
theory  for  laminated  beams  in  which  the  equilibrium  equations 
are  assumed  to  be  satisfied  in  an  average  sense  over  the  width 
of  the  beam.  A  new  set  of  force  and  moment  resultants  for  the 
beam  were  introduced  from  this  assumption.   Two  problems  of 
practical  interest  were  solved  using  this  method.  One  of  them 
was    the    problem    of    a    specially    orthotropic  delaminated 
cantilever  beam  under  torsion  loading,  and  the  other  was  the 
case  of  specially  orthotropic  and  quasi-orthotropic  simply 
supported  beams  under  transverse  loading.  In  both  problems  the 
delamination  was  assumed  to  be  in  different  positions  along 
the  thickness  of  the  beam,  and  the  strain  energy  release  rate 
for  each  position  was  computed.  The  strain  energy  release  rate 
was  found  to  have  the  maximum  value  in  the  middle  plane  of  the 
beam,    for  both  the  cantilever   beam  under   torsion   and  the 
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simply  supported  beam  under  transverse  force.  Practically, 
this  indicates  that  the  probability  of  crack  propagation 
becomes  higher  as  the  position  of  the  delamination  approaches 
the  middle  plane  of  the  beam. 

5.2  Future  research 

The  present  method  has  almost  all  its  results  validated 
by  comparison  with  other  methods  and  approaches.  The  method 
needs  to  be  applied  in  an  extensive  body  of  experimental 
tests. 

The  extension  of  this  method  to  the  study  of  buckling  in 
laminated  composite  beams  will  add  an  important  possibility  in 
composite  structural  mechanics. 

One  suggestion  that  naturally  arises  from  this  study  is 
to  extend  the  presently  developed  method  to  three-dimensional 
fracture  problems,  to  compute  strain  energy  release  rate  in 
delaminated  plates.  The  study  of  plates  with  multiple 
delaminations  can  be  an  important  sequel  to  this  research. 

The  final  suggestion  is  the  extension  of  this  method  to 
dynamic  problems,  which  will  bring  significant  benefits  to 
both  research  in  impact  dynamics  and  design  of  composite 
materials.  Then,  this  work  will  have  achieved  an  important 
goal  in  the  process  of  simulation  of  delaminations  and  crack 
propagation. 


APPENDIX 
STIFFNESS  MATRIX 


The  subroutine  STIFF  is  used  to  form  the  stiffness 
matrix  K  of  the  beam  finite  element  used  in  modeling 
anisotropic  laminated  beams.  The  matrix  K  results  from  the 
addition  of  two  matrices,  K  and  K. 

K  -  K  +  K 

Figure  A  illustrates  the  formation  process  of  matrix  K 
and  it  is  important  for  understanding  the  subroutine  steps. 
The  final  result  is  the  matrix  C,  with  dimension  of  14  x  14. 
The  input  to  the  subroutine  STIFF  are: 
B:  width  of  the  beam 
SV:  vector  with  the  stiffness  terms  for  top,  bottom  and 
regular  elements. 
XLE:  element  length 

IXl:  integer  that  defines  the  element  type:  l  for  top 

element,  2  for  bottom  element  and  3  for  regular  element. 
The  output  of  the  subroutine  STIFF  are: 
CC:  matrix  (7,7)    (see  Fig.  A.l) 
CA:  matrix  (7,14) (see  Fig.  A.l) 
C:  matrix  (14,14),  representative  of  the  stiffness  of  the 
beam  finite  element. 
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SUBROUTINE  STIFF 
Input  data 


matrix  k 

CHAT') 


14 
15 


21 


1                             14  15  21 

CB 

ACB 

CAB 

CCB 

1                               14  15  21 

CH 

ACH 

CAH 

CCH 

1                              14  15  21 

C 

AC 

CA 

CC 

120 


14 


14 


Figure  A  Stiffness 


matrix  K 


SUBROUTINE  STIFF ( B , SV , XLE , IXl , CC , CA , C ) 

IMPLICIT  DOUBLE  PRECISION (A-H, 0-Z) 

DIMENSION  SV(30) ,CB (14,14) ,CH( 14,14) ,C (14, 14) 

DIMENSION  CI (14, 7) ,C2 (14, 14) ,CC(7,7) , AC (14, 7) ,CA(7,14) 

DIMENSION  ACB(14,7) ,ACH(14,7) , CAB (7, 14) ,CAH(7,14) 

DIMENSION  CCB(7,7) ,CCH(7,7) 

11=10* (IXl-1) 

A11=SV(II+1) 

A16=SV(II+2) 

A66=SV(II+3) 

Bll=SV(II+4) 

B16=SV(II+5) 

B66=SV(II+6) 

Dll=SV(II+7) 

D16=SV(II+8) 

D66=SV(II+9) 

A55=SV(II+10) 

DO  1  1=1,14 

DO  1  J=l,14 

C(I, J)=0.0 

CONTINUE 

DO  5  1=1,14 

DO  5  J=l,14 

CB(I, J)=0.0 

CH(I, J)=0.0 

C2(I,J)=0.0 

CONTINUE 

DO  10  1=1,7 

DO  10  J=l,7 

CC(I,J)=0.0 

CCB(I,J)=0.0 

CCH(I, J)=0.0 

CONTINUE 

DO  15  1=1,14 

DO  15  J=l,7 

AC(I,J)=0.0 

CA(J,I)=0.0 

ACB(I, J)=0.0 

CAB (J, I) =0.0 

ACH(I, J)=0.0 

CAH(J,I)=0.0 

C1(I, J)=0.0 

CONTINUE 

'BAR'  MATRICES  

CB(14,14)  

CB(1,1)=A11*(70./XLE) 

CB(1,2)=A16*(70./XLE) 

CB(1,4)=A16*(-15.) 

CB(1,5)=B11*(70./XLE) 

CB(1,6)=B16*(-15. ) 

CB(l,7)=B16*(-70./XLE) 

CB(1,8)=A11*10./XLE 


CB(1,9)=A16*10./XLE 

CB(1,11)=A16*5. 

CB(1,12)=B11*10. /XLE 

CB(1,13)=B16*5. 

CB(1,14)=B16*(-10./XLE) 

CB(2,2)=A66*(70./XLE) 

CB(2,4)=A66*(-15. ) 

CB(2,5)=-CB(1,7) 

CB(2,6)=B66*(-15) 

CB(2,7)=B66*(-70./XLE) 

CB(2,8)=CB(1,9) 

CB(2,9)=A66*10./XLE 

CB(2, 11)=A66*5. 

CB(2,12)=-CB(1,14) 

CB(2, 13)=B66*5. 

CB(2,14)=B66*(-10./XLE) 

CB(3,3)=A55*(70./XLE) 

CB(3,5)=A55*(-15.) 

CB ( 3 , 10) =A55*10 . /XLE 

CB(3, 12)=A55*5. 

CB(4,4)=A66*4.*XLE 

CB(4,5)=CB(1,6) 

CB(4,6)=B66*4.*XLE 

CB(4,7)=-CB(2,6) 

CB(4,8)=-CB(1,11) 

CB(4,9)=-CB(2,11) 

CB(4, 11)=A66*(-XLE) 

CB(4,12)=-CB(1,13) 

CB (4 , 13 ) =B66* (-XLE) 

CB(4,14)=CB(2,13) 

CB(5, 6)=D16*(-15) 

CB (5, 5) =D11* (70 . /XLE) +A55*4 . *XLE 

CB(5,7)=D16*(-70./XLE) 

CB(5,8)=CB(1,12) 

CB(5,9)=-CB(1,14) 

CB(5,10)=-CB(3,12) 

CB(5,11)=CB(1,13) 

CB(5,12)=D11*10./ XLE+A5  5  * ( -XLE ) 

CB(5,13)=D16*5. 

CB(5,14)=D16*(-10./XLE) 

CB(6,6)=D66*4.*XLE 

CB(6,7)=D66*15 

CB(6,8)=-CB(1,13) 

CB(6,9)=-CB(2,13) 

CB(6,11)=CB(4,13) 

CB(6,12)=-CB(5,13) 

CB(6,13)=D66*(-XLE) 

CB(6,14)=D66*5. 

CB(7,7)=D66*70./XLE 

CB(7,8)=CB(1,14) 

CB(7,9)=CB(2,14) 

CB(7,11)=-CB(2,13) 


CB(7,12)=CB(5,14) 

CB(7,13)=-CB(6,14) 

CB (7 , 14) =D66*10 . /XLE 

CB(8,8)=A11*70. /XLE 

CB(8,9)=A16*70. /XLE 

CB(8, 11)=A16*15. 

CB (8 , 12 ) =B11*70 . /XLE 

CB(8,13)=B16*15. 

CB(8, 14)=CB(1,7) 

CB(9,9)=CB(2,2) 

CB(9,11)=-CB(2,4) 

CB(9,12)=-CB(1,7) 

CB(9,13)=-CB(2,6) 

CB(9,14)=CB(2,7) 

CB(10,10)=CB(3,3) 

CB(10,12)=-CB(3,5) 

CB(11,11)=CB(4,4) 

CB(11,12)=-CB(1,6) 

CB(11, 13)=CB(4,6) 

CB(11, 14)=CB(2,6) 

CB ( 12 , 12 ) =D11*70 . /XLE+A55*4 . *XLE 

CB(12, 13)=D16*15, 

CB{12,14)=CB(5,7) 

CB(13,13)=CB(6,6) 

CB(13,14)=-CB(6,7) 

CB(14,14)=CB(7,7) 

DO  20  1=1,14 

DO  20  J=I,14 

CB(J,I)=CB(I, J) 

CONTINUE 

DO  18  1=1,14 

DO  18  J=l,14 

CB(I, J)=CB(I, J) *B/30. 

CONTINUE 

MATRIX  ACB(14,7)  


ACB(1 

,  1 

)=All*(-80. 

/XLE) 

ACB(1 

r2 

l=A16*(-80. 

/XLE) 

ACB(1 

,4 

i=A16*(-20. 

) 

ACB{1 

r5 

=Bll*(-80. 

/XLE) 

ACB(1 

,6 

)=B16*(-20. 

) 

ACB(1 

=B16*(80./XLE) 

ACB(2, 

i; 

=ACB (1,2) 

ACB(2, 

2 

=A66*(-80. 

/XLE) 

ACB(2, 

4, 

=A66*(-20. 

) 

ACB(2, 

5] 

=-ACB (1,7) 

ACB(2, 

6] 

=B66*(-20. 

) 

ACB ( 2 , 

7] 

=B66*(80./XLE) 

ACB ( 3 , 

3] 

=A55*(-80. 

/XLE) 

ACB ( 3 , 

5) 

=A55*(-20. 

) 

ACB (4, 

1) 

=A16*20. 

ACB(4,2)=A66*20. 
ACB (4 , 4) =A66*2 . *XLE 
ACB(4,5)=B16*20. 
ACB(4 , 6) =B66*2 . *XLE 
ACB(4,7)=B66*(-20.) 
ACB{5,1)=ACB(1,5) 
ACB(5,2)=-ACB(1,7) 
ACB(5, 3)=A55*20. 
ACB(5,4)=-20.*B16 

ACB (5,5) =D11* (-80 . /XLE) +A55*2 . *XLE 

ACB(5,6)=D16*(-20,) 

ACB(5,7)=D16*80. /XLE 

ACB(6, 1)=ACB(4,5) 

ACB(6,2)=-ACB(4,7) 

ACB(6,4)=ACB(4,6) 

ACB(6,5)=D16*20. 

ACB (6 , 6) =D66*2 . *XLE 

ACB(6,7)=D66*(-20. ) 

ACB(7,1)=ACB(1,7) 

ACB(7,2)=ACB(2,7) 

ACB(7,4)=-ACB(2,6) 

ACB(7,5)=ACB(5,7) 

ACB(7,6)=D66*(20. ) 

ACB (7 , 7 ) =D66* (-80 . /XLE) 

ACB(8,1)=ACB(1,1) 

ACB(8,2)=ACB(1,2) 

ACB(8,4)=ACB(4,1) 

ACB(8,5)=ACB(1,5) 

ACB(8,6)=ACB(4,5) 

ACB(8,7)=ACB(1,7) 

ACB(9,1)=ACB(1,2) 

ACB(9,2)=ACB(2,2) 

ACB(9,4)=ACB(4,2) 

ACB(9,5)=-ACB(1,7) 

ACB(9,6)=-ACB(4,7) 

ACB(9,7)=ACB(2,7) 

ACB(10,3)=ACB(3,3) 

ACB(10,5)=A55*20. 

ACB(11,1)=-ACB(4,1) 

ACB(11,2)=-ACB(4,2) 

ACB(11,4)=ACB(4,4) 

ACB(11,5)=-ACB(4,5) 

ACB(11,6)=ACB(4,6) 

ACB(11,7)=20.*B66 

ACB(12,1)=ACB(1,5) 

ACB(12,2)=-ACB(1,7) 

ACB(12,3)=-ACB(5,3) 

ACB(12,4)=20.*B16 

ACB ( 12 , 5) =D11* (-80 . /XLE) +A55*2 . *XLE 

ACB(12,6)=ACB(6,5) 

ACB(12,7)=ACB(5,7) 


ACB(13, 1)=-ACB(4,5) 

ACB(13,2)=ACB(4,7) 

ACB(13,4)=ACB(4,6) 

ACB(13,5)=-ACB(6,5) 

ACB (13,6) =D66*2 . *XLE 

ACB(13,7)=D66*(20.) 

ACB(14,1)=ACB(1,7) 

ACB(14,2)=ACB(2,7) 

ACB(14,4)=ACB(4,7) 

ACB(14,5)=ACB(5,7) 

ACB(14,6)=ACB(6,7) 

ACB(14,7)=ACB(7,7) 

DO  24  1=1,14 

DO  24  J=l,7 

ACB(I, J)=ACB(I, J) *B/30. 

CONTINUE 
.CAB (7, 14)  

DO  25  1=1,14 

DO  25  J=l,7 

CAB(J,I)=ACB(I,J) 

CONTINUE 
.CC(7,7)  

CCB(1, 1)=A11*160.  /XLE 

CCB (1,2) =A16*160 . /XLE 

CCB(1,5)=B11*160./XLE 

CCB(1,7)=-B16*160./XLE 

CCB(2,2)=A66*160. /XLE 

CCB(2,5)=-CCB(1,7) 

CCB(2,7)=B66*(-160./XLE) 

CCB(3,3)=A55*160. /XLE 

CCB(4,4)=A66*16.*XLE 

CCB(4, 6)=B66*16. *XLE 

CCB (5 , 5) =011*160 . /XLE+A55*16 . *XLE 

CCB(5,7)=D16*(-160./XLE) 

CCB(6, 6)=D66*16.*XLE 

CCB (7 , 7) =066*160 . /XLE 

DO  28  1=1,7 

DO  28  J=I,7 

CCB(I, J)=CCB(I,J)*B/30. 

CONTINUE 

DO  26  1=1,7 

DO  26  J=I,7 

CCB(J,I)=CCB(I,  J) 

CONTINUE 

'HAT'  MATRICES   

CH(14,14)  

BB=B**3/12. 

CH(4,4)=A11*70./XLE 

CH(4,6)=B11*70./XLE 

CH(4,11)=A11*10./XLE 

CH(4,13)=B11*10./XLE 


CH(6,6)=(D11*70./XLE)+(A55*4.*XLE) 

CH(6,7)=A55*(-15.) 

CH(6,11)=CH(4,13) 

CH(6,13)=(D11*10./XLE)+A55*(-XLE) 
CH(6,14)=A55*{-5.) 

CH(7,7)=A55*70./XLE 
CH(7,13)=-CH(6,14) 
CH(7,14)=A55*10./XLE 
CH(11,11)=CH(4,4) 
CH(11,13)=CH(4,6) 
CH(13,13)=CH(6,6) 
CH(13,14)=-CH(6,7) 
CH(14,14)=CH(7,7) 
DO  55  1=1,14 
DO  55  J=I,14 
CH(J,I)=CH(I,J) 
CONTINUE 
DO  60  1=1,14 
DO  60  J=l,14 
CH(I, J)=CH(I, J) *BB/30. 
CONTINUE 
.ACH(14,7)  

ACH(4,4)=BB*(All*(-80./XLE))/30. 

ACH(4,6)=BB*(Bll*(-80./XLE))/30. 

ACH(6,4)=BB*(Bll*(-80./XLE) ) /30. 

ACH(6 , 6) =BB* (Dll* (-80 . ) /XLE+A55*2 . *XLE) /30 

ACH(6,7)=BB*(A55*20.)/30.  ^^^/-^u. 

ACH(7,6)=BB*{A55*(-20.) ) /30. 

ACH(7,7)=BB*(A55*(-80./XLE))/30. 

ACH{ll,4)=BB*(All*(-80./XLE))/30. 

ACH(ll,6)=BB*(Bll*(-80./XLE))/30. 

ACH(13,4)=BB*(Bll*(-80./XLE))/30. 

ACH(13 , 6) =BB* (Dll* (-80 . ) /XLE+A55*2 . *XLE) /30 

ACH(13,7)=-BB*(A55*20.)/30.  ' 

ACH(14,6)=BB*(A55*20.) /30. 

ACH(14,7)=BB*(A55*(-80./XLE) ) 730. 
.CAH(7,14)  ' 

DO  65  1=1,14 

DO  65  J=l,7 

CAH(J,I)=ACH(I, J) 

CONTINUE 

.CCH(7,7)  

CCH(4,4)=BB*(A11*160.'/XLE) /30. 
CCH(4,6)=BB*(B11*160./XLE) /30. 

CCH(6,6)=BB*( (011*160. /XLE)+(A55*16.*XLE) ) 730 
CCH(7,7)=BB*(A55*160./XLE)/30.  ^LE))/30, 


DO  70  1=1,7  127 
DO  70  J=I,7 
CCH(J,I)=CCH(I, J) 
70  CONTINUE 

C  ADDING  THE   'BAR'  AND   'HAT'  MATRICES  

DO  30  1=1,14 

DO  30  J=l,14 

C(I, J)=CB(I, J)+CH(I, J) 
3  0  CONTINUE 

DO  32  1=1,14 

DO  32  J=l,7 

AC(I, J)=ACB(I, J)+ACH(I, J) 
CA(J,I)=CAB(J,I)+CAH(J,I) 
32  CONTINUE 

DO  34  1=1,7 
DO  34  J=l,7 

CC(I, J)=CCB(I, J)+CCH(I, J) 
34  CONTINUE 

C  STATIC  CONDENSATION  

CALL  INVDET(CC,7) 

CALL  MATMUL(AC,CC,C1,14,7,7) 

CALL  MATMUL(C1,CA,C2,14,7,14) 

DO  40  1=1,14 

DO  40  J=l,14 

C(I,J)=C(I,J)-C2(I,J) 
40  CONTINUE 

CALL  INVDET(CC,7) 

RETURN 

END  - 
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